H2O CVRQD potential (#171)

Single Water PES published in:
"CVRQD ab initio ground-state adiabatic
potential energy surfaces for the water molecule"

https://doi.org/10.1063/1.2378766

.fprettify.rc

@@ -26,4 +26,4 @@ enable-decl=True
 disable-fypp=True
 
 case=[1,1,1,2]
-exclude=[random.F90,fftw3.F90,force_cp2k.F90, h2o_schwenke.f]
+exclude=[random.F90,fftw3.F90,force_cp2k.F90, h2o_schwenke.f, h2o_cvrqd.f]

src/Makefile

@@ -1,6 +1,6 @@
 F_OBJS := constants.o fortran_interfaces.o error.o modules.o mpi_wrapper.o files.o utils.o io.o random.o arrays.o qmmm.o fftw_interface.o \
           shake.o nosehoover.o gle.o transform.o potentials.o force_spline.o estimators.o ekin.o vinit.o plumed.o \
-          remd.o force_bound.o water.o h2o_schwenke.o force_h2o.o force_cp2k.o sh_integ.o surfacehop.o landau_zener.o\
+          remd.o force_bound.o water.o h2o_schwenke.o h2o_cvrqd.o force_h2o.o force_cp2k.o sh_integ.o surfacehop.o landau_zener.o\
           force_mm.o tera_mpi_api.o force_abin.o force_tcpb.o force_tera.o force_terash.o en_restraint.o analyze_ext_template.o geom_analysis.o analysis.o \
           minimizer.o mdstep.o forces.o cmdline.o init.o
 

src/force_h2o.F90

@@ -57,6 +57,9 @@ subroutine force_h2o(x, y, z, fx, fy, fz, eclas, natom, nbeads)
       ! TODO: Use function pointers to select the potential and pass it down
       if (h2opot == 'schwenke') then
          call force_h2o_schwenke(x, y, z, fx, fy, fz, eclas, natom, nbeads)
+      else if (h2opot == 'cvrqd') then
+         call force_h2o_cvrqd(x, y, z, fx, fy, fz, eclas, natom, nbeads)
+         call fatal_error(__FILE__, __LINE__, 'Numerical forces not yet implemented!')
       else
          call fatal_error(__FILE__, __LINE__, 'Potential '//trim(h2opot)//' not implemented')
       end if
@@ -97,6 +100,39 @@ subroutine force_h2o_schwenke(x, y, z, fx, fy, fz, Eclas, natom, nbeads)
 
    end subroutine force_h2o_schwenke
 
+   subroutine force_h2o_cvrqd(x, y, z, fx, fy, fz, Eclas, natom, nbeads)
+      real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
+      real(DP), intent(inout) :: fx(:, :), fy(:, :), fz(:, :)
+      real(DP), intent(inout) :: Eclas
+      integer, intent(in) :: natom, nbeads
+      ! Internal water coordinates
+      real(DP) :: rOH1, rOH2, aHOH_rad
+      real(DP) :: E
+      integer :: iw
+      real(DP) :: mH, mO
+
+      ! TODO: Pass in the actual masses
+      ! What should be these set to? Should we set pure isotopes?
+      mH = 1.008D0
+      mO = 15.999D0
+
+      Eclas = 0.0D0
+      ! The H2O potentials are evaluated in internal coordinates, but ABIN works in cartesians
+      do iw = 1, nbeads
+         call get_internal_coords(x, y, z, iw, rOH1, rOH2, aHOH_rad)
+
+         call h2o_pot_cvrqd(E, rOH1, rOH2, aHOH_rad, mO, mH)
+
+         Eclas = Eclas + E
+      end do
+      Eclas = Eclas / nbeads
+
+      ! TODO: Given the small difference between the Schwenke potential,
+      ! we might not need to implement numerical forces here.
+      ! call numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads)
+
+   end subroutine force_h2o_cvrqd
+
    ! TODO: Implement numerical forces generally for all potentials
    ! For now, they can be implemented here and hardcoded for a specific H2O potential
    subroutine numerical_forces(x, y, z, Epot, fx, fy, fz, natom, nbeads)

src/fortran_interfaces.F90

@@ -63,14 +63,21 @@ function usleep(useconds) bind(c, name='usleep')
          integer(kind=C_INT) :: usleep
       end function usleep
 
-      ! Returns a potential energy of a water molecule
+      ! Computes potential energy of a water molecule
       ! using Schwenke potential, see h2o_schwenke.f
       subroutine h2o_pot_schwenke(rij, v, n)
          import :: DP
          integer, intent(in) :: n
-         real(DP) :: rij(n, 3), v(n)
+         real(DP), intent(in) :: rij(n, 3)
+         real(DP), intent(out) :: v(n)
       end subroutine h2o_pot_schwenke
 
+      subroutine h2o_pot_cvrqd(V, rOH1, rOH2, aHOH, mH, mO)
+         import :: DP
+         real(DP), intent(out) :: V
+         real(DP), intent(in) :: rOH1, rOH2, aHOH, mH, mO
+      end subroutine h2o_pot_cvrqd
+
    end interface
 
 end module mod_interfaces