Code design for QMD simulations

Hi everyone,

Thank you in advance for taking the time to read my post. I apologise for bringing C code to a Fortran community :dagger: but it is necessary to clarify question.

I am currently re-writing my C code for running Quantum Molecular Dynamics simulations in Fortran. In my C code I have essentially taken an OOP approach whereby I define a header file potential.h

#ifndef POTENTIAL_DOT_H
#define POTENTIAL_DOT_H

# include "complex.h"

struct Potential;
struct Potential{
  double delta; // parameters
  double alpha; 
  double eps;
  double gamma;
  double (*func_v)(struct Potential *pot, double *x, unsigned int dim);
  void (*func_v12d)(struct Potential *pot, double *x, double *grad, unsigned int dim);
  double complex (*func_get_tau)(struct Potential *pot);
};

/* -------------- Construct potential from model systems -------- */

struct Potential    *potential_construct(
    double (*func_v)(struct Potential *pot, double *x, unsigned int dim),
    void (*func_v12)(struct Potential *pot, double *x, double *grad, unsigned int dim),
    double complex (*func_get_tau)(struct Potential *pot),
    char* potential_name, double param[]);

/* functions to be defined in the model potential*/
double v_v(struct Potential *pot, double *x, unsigned int dim);
void v_vd(struct Potential *pot, double *x, double *grad, unsigned int dim);
/* */
double complex get_tau(struct Potential *pot);
#endif

And then anytime I want to define a new potential I simply include the
header file “potential.h”. For instance, for a Landau-Zener potential:

#include <math.h>
#include "potential.h"

/*
* The Landau-Zener type potential is given in ?
*/

/* Handles construction of potential in the adiabatic representation */

double v_v(struct Potential *pot, double *x, unsigned int dim){
 return pot->delta;
}

double v_vd(struct Potential *pot, double *x, unsigned int dim){
 return 0.0;
}

double get_tau(struct Potential *pot){
 return M_PI * pow(pot->delta,2) / 2 / pot->alpha;
}

I essentially want to do the same in my Fortran code. But which approach is best? An OOP approach or something else? And how should I design the code?

3 Likes

@MichaelRedenti , you will receive diverse opinions here that will help you decide on your approach.

But very quickly, as things stand with Fortran currently where every derived type (struct / class) is an extensible type, you will find the OOP approach to be helpful

  1. if inheritance and polymorphism are important in your overall simulation framework,
  2. with the notion of a namespace / packaging or an assembly since the language itself, like C, does not include such a facility intrinsically,
  3. to have a near one-to-one correspondence with your current code in C, at least conceptually.

If you instead went with a Fortran module-oriented object-based approach, you may gain some in performance on account of avoiding dynamic dispatch at run-time but at a slight cost of having to generally do object composition vs inheritance.

Either way, you are asking the right questions; moreover, you are absolutely on the right track in considering Fortran for your QMD simulations. You will find a lot of discussion around marketing of Fortran at this Discourse - in your case, if I may paraphrase a TV commercial - “you will like the way your QMD simulations code is gonna look in Fortran, I guarantee it”!!

2 Likes

@FortranFan , would the Fortran module-oriented object-based approach consist of modules and submodules. For instance, I could define a module potential

!> Potential base module
!>
!> ...
module potential
  implicit none
  private
 
public :: pot

  !> Evaluate potential at point x
  interface pot
    pure module function pot_v(x) result(v)
      real(kind), intent(in) :: x
      real(dp) :: v
    end function pot_v
    module function pot_vd(x, dim) result(vd)
      real(kind), intent(in) :: x(:)
      integer, intent(in) :: dim
      real(dp), dimension(dim) :: vd
    end function pot_vd
  end interface pot

end module potential

and then an actual potential would be implemented through a submodule as

submodule (potential) potential_lz
  implicit none

contains

  pure module function pot_v(x) result(v)
    real(kind), intent(in) :: x
    real(kind) :: v

   !... 
  end function pot_v

! similarly for pot_vd
  ! ...
end submodule potential_lz
1 Like

Yes, that’s one way to about it. You can define a derived type also to go with and have dummy arguments of the same type with your potential functions.

   ..
   type :: Potential_t
      real(kind=..) :: delta = <default value> ! parameters
      real(kind=..) :: alpha = <default value>
      real(kind=..) :: eps = <default value>
      real(kind=..) :: gamma = <default value>
   ..
   end type
   ..
   interface pot
      pure module function pot_v(pot, x) result(v)
         type(Potential_t), intent(in) :: pot
1 Like

Note if this is a common situation, then it implies inheritance and polymorphism might be important for you and in that case, your initial thought of OOP approach can work well.

1 Like

Although a Landau-Zener potential (and any other potential) would be an instantiation of the potential class rather than a “child”. So this would not really count as inheritance would it?

Well, you will know Fortran allows you to define an ABSTRACT type (type, abstract), say for your potential “class” that has the base fields for such a class and with virtual methods.

You can then have a Landau-Zener potential type that has a concrete definition of such an abstract class with actual methods toward Landau-Zener potential functions that override the virtual ones.

This OOP approach seems to be a better fit for what I am trying to do compared to the module-submodule way (I am not sure about the performance differences or what happens under the hood).

I am saying this simply because if I were to take the module-submodule approach then I would have to pass both type(Potential_t) :: params and the pot_v(pot,x) function to say a numerical solver each time. Whereas with the OOP approach attribute and methods would simply be part of the class object and what would have to pass only one thing making life easier.

I’m developing in a library of potential for quantum dynamic in the diabatic or adiabatic representation (see the link below). It uses OOP and Automatic differentiation to get the gradient and the hessian, so that implementation of potential is fast:

  • a modification of the initialization procedure.
  • a file with the potential expression
1 Like

@MichaelRedenti , that’s what I was referring to in my first response upthread, see point 2 which includes the benefits you notice.

You can make things work either way, OOP is not a must. But you’re right indeed that for those who are so inclined, OOP can offer particular advantages with structured programming that can make certain aspects easier. You’ll find though that is not the case with everyone else - there are others who find the OOP approach itself tedious and therefore eschew it, as you will know. Hence the need for a balanced view, particularly if you intend to collaborate with others on this code for QMD simulations.

1 Like