Integrate the Schwenke potential

.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]
+exclude=[random.F90,fftw3.F90,force_cp2k.F90, h2o_schwenke.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 force_cp2k.o sh_integ.o surfacehop.o landau_zener.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\
           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
 
@@ -27,5 +27,8 @@ clean:
 %.o: %.F90
 	$(FC) $(FFLAGS) $(DFLAGS) -c $<
 
+%.o: %.f
+	$(FC) $(FFLAGS) $(DFLAGS) -c $<
+
 %.o: %.cpp
 	$(CXX) $(CXXFLAGS) $(DFLAGS) -c $<

src/force_abin.F90

@@ -5,7 +5,7 @@
 ! The basic workflow is very simple:
 ! 1. ABIN writes current geometry into file geom.dat
 !   in a XYZ format without the header.
-! 2. ABIN launches the shell script POT/r.pot
+! 2. ABIN launches the shell script <POT>/r.<pot>
 ! 3. The shellscript does a few things:
 !    i) Takes the input geometry and prepares the input file
 !    ii) Launches the QM program
@@ -15,7 +15,7 @@
 !
 ! NOTE: We append bead index to every file name so that we can
 ! call the interface in parallel in PIMD simulations.
-! NOTE: Interface for Surface Hopping is a bit more complicated,
+! NOTE: The interface for Surface Hopping is a bit more complicated,
 ! see interfaces/MOLPRO-SH/r.molpro-sh
 module mod_shell_interface_private
    use mod_const, only: DP, ANG

src/force_h2o.F90

@@ -0,0 +1,53 @@
+! This is invoked via pot='_h2o_'
+module mod_h2o_pes
+   use mod_const, only: DP
+   use mod_error, only: fatal_error
+   implicit none
+   private
+   public :: force_h2o_schwenke
+   save
+contains
+
+   subroutine force_h2o_schwenke(x, y, z, fx, fy, fz, eclas, natom, walkmax)
+      use mod_utils, only: get_distance, get_angle
+      use mod_const, only: PI
+      real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
+      real(DP), intent(inout) :: fx(:, :), fy(:, :), fz(:, :)
+      real(DP), intent(inout) :: eclas
+      integer, intent(in) :: natom, walkmax
+      ! Internal water coordinates
+      real(DP) :: rOH1, rOH2, aHOH_deg, aHOH_rad
+      real(DP) :: rij(walkmax, 3)
+      real(DP) :: Epot(walkmax)
+      integer :: iw
+
+      ! TODO: Check atom names
+      if (natom /= 3) then
+         call fatal_error(__FILE__, __LINE__, "This is not a water molecule!")
+      end if
+
+      ! Compute internal coordinates from cartesians
+      ! OH distances in bohrs, HOH angle in radians
+      do iw = 1, walkmax
+         rOH1 = get_distance(x, y, z, 1, 2, iw)
+         rOH2 = get_distance(x, y, z, 1, 3, iw)
+         aHOH_deg = get_angle(x, y, z, 2, 1, 3, iw)
+         aHOH_rad = aHOH_deg * PI / 180.0D0
+         print*, 'Angle = ', aHOH_deg, aHOH_rad
+         rij(iw, 1) = rOH1
+         rij(iw, 2) = rOH2
+         rij(iw, 3) = aHOH_rad
+      end do
+
+      call pes_h2o_schwenke(rij, Epot, walkmax)
+
+      ! TODO: Implement numerical forces
+      fx = 0.0D0; fy = 0.0D0; fz = 0.0D0
+
+      do iw = 1, walkmax
+         eclas = eclas + Epot(iw)
+      end do
+      eclas = eclas / walkmax
+   end subroutine force_h2o_schwenke
+
+end module mod_h2o_pes

src/forces.F90

@@ -128,6 +128,7 @@ subroutine force_wrapper(x, y, z, fx, fy, fz, e_pot, chpot, walkmax)
    use mod_const, only: DP
    use mod_general, only: natom, ipimd
    use mod_water, only: watpot
+   use mod_h2o_pes, only: force_h2o_schwenke
    use mod_force_mm, only: force_mm
    use mod_sbc, only: force_sbc, isbc
    use mod_plumed, only: iplumed, force_plumed
@@ -155,6 +156,8 @@ subroutine force_wrapper(x, y, z, fx, fy, fz, e_pot, chpot, walkmax)
       call force_mm(x, y, z, fx, fy, fz, eclas, walkmax)
    case ("_mmwater_")
       call force_water(x, y, z, fx, fy, fz, eclas, natom, walkmax, watpot)
+   case ("_h2o_")
+      call force_h2o_schwenke(x, y, z, fx, fy, fz, eclas, natom, walkmax)
    case ("_splined_grid_")
       ! Only 1D spline grid supported at the moment
       call force_splined_grid(x, fx, eclas, walkmax)

src/fortran_interfaces.F90

@@ -63,6 +63,14 @@ function usleep(useconds) bind(c, name='usleep')
          integer(kind=C_INT) :: usleep
       end function usleep
 
+      ! Returns a potential energy of a water molecule
+      ! using Schwenke potential, see h2o_schwenke.f
+      subroutine pes_h2o_schwenke(rij, v, n)
+         import :: DP
+         integer, intent(in) :: n
+         real(DP) :: rij(n, 3), v(n)
+      end subroutine pes_h2o_schwenke
+
    end interface
 
 end module mod_interfaces