Add Schwenke's gas phase H2O analytical potential (#166)

.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]

.github/workflows/gfortran.yml

@@ -29,18 +29,17 @@ env:
 jobs:
 
   format:
-    name: Format check
-    runs-on: ubuntu-20.04
-    strategy:
-      fail-fast: false
+    name: Code format check
+    runs-on: ubuntu-latest
 
     steps:
     - uses: actions/checkout@v4
+    - uses: actions/setup-python@v5
       with:
-        fetch-depth: 2
+        python-version: '3.11'
 
     - name: Install fprettify dependencies
-      run: pip install --user configargparse
+      run: pip install configargparse
 
     - name: Run fprettify
       run: |

autoformat.py

@@ -31,8 +31,8 @@
 try:
     import configargparse
 except ImportError:
-    print("Could not find configargparse package")
-    print("Please install it e.g. via 'pip install --user configargparse")
+    print("ERROR: Could not find configargparse package")
+    print("Please install the package, e.g. 'pip install --user configargparse'")
     sys.exit(1)
 
 if len(sys.argv) == 1:

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,117 @@
+! Interface to analytical H2O (single water molecule) potentials.
+! This is invoked via pot='_h2o_', and the potential is selected via
+! h2opot='schwenke' in namelist General in ABIN input file.
+module mod_force_h2o
+   use mod_const, only: DP
+   use mod_files, only: stdout, stderr
+   use mod_error, only: fatal_error
+   use mod_interfaces, only: h2o_pot_schwenke
+   implicit none
+   private
+   public :: initialize_h2o_pot, force_h2o, h2opot
+   character(len=20) :: h2opot = 'schwenke'
+   save
+contains
+
+   ! Perform some basic validation
+   subroutine initialize_h2o_pot(natom, atom_names)
+      integer, intent(in) :: natom
+      character(len=2), intent(in) :: atom_names(natom)
+
+      if (natom /= 3) then
+         call fatal_error(__FILE__, __LINE__, "This is not a water molecule!")
+      end if
+
+      if (atom_names(1) /= 'O' .or. atom_names(2) /= 'H' .or. atom_names(3) /= 'H') then
+         write (stderr, *) 'ERROR: Bad element type.'
+         write (stderr, *) 'Water atoms must be ordered as "O H H"'
+         call fatal_error(__FILE__, __LINE__, "This is not a water molecule!")
+      end if
+      ! TODO: Some PES initialization could be performed here
+   end subroutine initialize_h2o_pot
+
+   ! Compute internal coordinates from cartesians,
+   ! OH distances in Bohrs, HOH angle in radians.
+   subroutine get_internal_coords(x, y, z, iw, rOH1, rOH2, aHOH_rad)
+      use mod_const, only: PI
+      use mod_utils, only: get_distance, get_angle
+      real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
+      integer, intent(in) :: iw
+      real(DP), intent(out) :: rOH1, rOH2, aHOH_rad
+      real(DP) :: aHOH_deg
+
+      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
+   end subroutine get_internal_coords
+
+   ! This is just a wrapper function that selects the H2O potential
+   ! based on user input. For now only h2opot="schwenke" is supported
+   subroutine force_h2o(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
+
+      ! 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
+         call fatal_error(__FILE__, __LINE__, 'Potential '//trim(h2opot)//' not implemented')
+      end if
+   end subroutine force_h2o
+
+   subroutine force_h2o_schwenke(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) :: rij(nbeads, 3)
+      real(DP) :: Epot(nbeads)
+      integer :: iw
+
+      ! 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)
+         rij(iw, 1) = rOH1
+         rij(iw, 2) = rOH2
+         rij(iw, 3) = aHOH_rad
+      end do
+
+      ! Calculated the potential for all PI beads at once
+      ! The potential is implemented in h2o_schwenke.f
+      call h2o_pot_schwenke(rij, Epot, nbeads)
+
+      ! For Path Integrals, the final energy of the PI necklace
+      ! is an average over all beads.
+      ! TODO: Forces need to be appropriately scaled as well!
+      do iw = 1, nbeads
+         Eclas = Eclas + Epot(iw)
+      end do
+      Eclas = Eclas / nbeads
+
+      call numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads)
+
+   end subroutine force_h2o_schwenke
+
+   ! 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)
+      real(DP), intent(in) :: x(natom, nbeads)
+      real(DP), intent(in) :: y(natom, nbeads)
+      real(DP), intent(in) :: z(natom, nbeads)
+      real(DP), intent(inout) :: fx(natom, nbeads)
+      real(DP), intent(inout) :: fy(natom, nbeads)
+      real(DP), intent(inout) :: fz(natom, nbeads)
+      ! This is the energy for the currrent geometry that has already been calculated
+      real(DP), intent(in) :: Epot(nbeads)
+      integer, intent(in) :: natom, nbeads
+
+      ! Need to implement numerical forces first
+      call fatal_error(__FILE__, __LINE__, 'Numerical forces not yet implemented!')
+   end subroutine numerical_forces
+
+end module mod_force_h2o

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_force_h2o, only: force_h2o
    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(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 h2o_pot_schwenke(rij, v, n)
+         import :: DP
+         integer, intent(in) :: n
+         real(DP) :: rij(n, 3), v(n)
+      end subroutine h2o_pot_schwenke
+
    end interface
 
 end module mod_interfaces