Schwenke numerical

.gitignore

@@ -10,3 +10,4 @@ make.vars
 bin/
 # from fprettify
 __pycache__/
+Analysis

sample_inputs/input.in.cmd

@@ -2,25 +2,25 @@ This is a sample input file for ABIN for Classical MD simulation --
 -- plain old Newton's equations of motion.
 
 &general
-pot='g09'		! where do we obtain forces and energies?
+pot='_h2o_'		! where do we obtain forces and energies?
 mdtype='md',		! md = classical MD
-nstep=50000,		! number of steps
+nstep=1000,		! number of steps
 dt=20.,			! timestep (a.u.)
 
 irest=0,		! should we restart from restart.xyz? (ignoring input coordinates and velocities)
 
 nwrite=1,		! how often some output should be printed (energies etc.)
                         ! 1 = every time step
                         ! 10 = each 10th time step
-nwritex=5,		! how often should we print coordinates?
+nwritex=1,		! how often should we print coordinates?
 nrest=5,		! how often we print restart files?
 nwritev=0,		! how often we print velocities? (by default ABIN does not print velocities)
 /
 
 
 &nhcopt
 temp=298.15,		! temperature [K] for Maxwell-Boltzmann sampling and thermostat
-inose=1,		! Thermostat options: 
+inose=0,		! Thermostat options: 
                         ! 0 = NVE( microcanonical)
                         ! 1 = Nose-Hoover
                         ! 2 = Generalized Langevin Equation (GLE)

src/force_h2o.F90

@@ -1,3 +1,4 @@
+
 ! 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.
@@ -99,19 +100,81 @@ 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)
+
+   subroutine numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads)
+      real(DP), intent(in) :: x(3, 1)
+      real(DP), intent(in) :: y(3, 1)
+      real(DP), intent(in) :: z(3, 1)
+      real(DP), intent(inout) :: fx(3, 1)
+      real(DP), intent(inout) :: fy(3, 1)
+      real(DP), intent(inout) :: fz(3, 1)
       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
+      ! Create new copies of arrays
+      real(DP) :: x_new_forward(3, 1)
+      real(DP) :: y_new_forward(3, 1)
+      real(DP) :: z_new_forward(3, 1)
 
-end module mod_force_h2o
+      ! This is the energy for the currrent geometry that has already been calculated
+      real(DP), intent(in) :: Epot(nbeads) !nbeads
+
+      ! Internal water coordinates
+      real(DP) :: new_rOH1, new_rOH2, new_aHOH_rad
+      real(DP) :: new_rij(1, 3)
+
+      ! Schwenke calculated peterbed geometry energy
+      real(DP) :: Epot_delta(nbeads) !nbeads
+
+      real(DP) :: Eclas_orig, Eclas_plus
+      real(DP) :: delta = 5.0E-5_DP
+      integer :: i, j, k, iw
+
+      ! Calculate forces numerically using central differences
+      do j = 1, nbeads !nbeads
+
+         ! Save the original energy
+         Eclas_orig = Epot(j)
+
+         do i = 1, natom !natom
+
+            do k = 1, 3 ! x, y, z
+
+               ! Copy the original atom coordinates
+               x_new_forward(:,:) = x(:,:)
+               y_new_forward(:,:) = y(:,:)
+               z_new_forward(:,:) = z(:,:)
+
+               ! Move the atom forwards
+               select case (k)
+                  case (1)
+                     x_new_forward(i,1) = x_new_forward(i,1) + delta
+                  case (2)
+                     y_new_forward(i,1) = y_new_forward(i,1) + delta
+                  case (3)
+                     z_new_forward(i,1) = z_new_forward(i,1) + delta
+               end select
+
+               ! Calculate the energy for the forward perturbed geometry
+               call get_internal_coords(x_new_forward, y_new_forward, z_new_forward, j, new_rOH1, new_rOH2, new_aHOH_rad)
+
+               new_rij(j, :) = [new_rOH1, new_rOH2, new_aHOH_rad]
+
+               call h2o_pot_schwenke(new_rij, Epot_delta, nbeads)
+
+               Eclas_plus = Epot_delta(j)
+
+               ! Calculate the numerical force
+               select case (k)
+                  case (1)
+                     fx(i,1) = -(Eclas_plus - Eclas_orig) / delta
+                  case (2)
+                     fy(i,1) = -(Eclas_plus - Eclas_orig) / delta
+                  case (3)
+                     fz(i,1) = -(Eclas_plus - Eclas_orig) / delta
+               end select
+            end do
+         end do
+      end do
+   end subroutine numerical_forces
+
+end module mod_force_h2o

src/h2o_scan.f90

@@ -0,0 +1,67 @@
+! To compile:
+! gfortran -c h2o_scan.f90
+! gfortran -fopenmp -L./ -L../water_potentials/ libabin.a ../water_potentials/libwater.a h2o_scan.o -labin -lwater -lm -lstdc++ -o scan
+program scan
+    !use mod_water
+    use mod_const, only: DP, ANG
+    use mod_init, only: read_xyz_file
+    implicit none
+    external :: force_water
+
+    integer, parameter :: NATOM = 3
+    integer, parameter :: NBEADS = 1
+    integer, parameter :: WATPOT = 1
+
+    ! Displacement for H atom at each timestep
+    real(DP) :: delta = 0.0005_DP
+    integer, parameter :: NSTEPS = 5000
+
+    real(DP) :: x(3, 1), y(3, 1), z(3, 1)
+    real(DP) :: fx(3, 1), fy(3, 1), fz(3, 1)
+    real(DP) :: energy, rOH
+    integer :: i, funit
+    character(len=2) :: atom_names(3)
+    character(len=256) :: fname
+
+    atom_names = (/ "O", "H", "H" /)
+
+    fname = "water.xyz"
+    open (newunit=funit, file=fname, action="read")
+    call read_xyz_file(funit, fname, atom_names, natom, nbeads, x, y, z)
+    close (funit)
+
+    ! Convert from Angstromgs to Bohr (atomic units)
+    x(:, 1) = x(:, 1) * ANG
+    y(:, 1) = y(:, 1) * ANG
+    z(:, 1) = z(:, 1) * ANG
+
+    ! print for output
+    print*, "Step   ", "rOH   ", "Energy"
+
+    open(99, file="water_geometries.xyz", action="write")
+    ! steps for stretching OH bond
+    do i = 1, NSTEPS
+        ! Calculate the initial OH bond length (rOH) from initial bond positions
+        rOH = sqrt((x(2, 1) - x(1, 1))**2 + (y(2, 1) - y(1, 1))**2 + (z(2, 1) - z(1, 1))**2)
+
+        ! Modify the OH bond length by changing the position of a single hydrogen atom
+        x(2, 1) = x(2, 1) + delta
+
+        ! Call the force_water function to calculate forces and energy
+        energy = 0.0D0
+        call force_water(x, y, z, fx, fy, fz, energy, natom, nbeads, watpot)
+
+        ! Print the current OH bond length and energy
+        print '(I4,F16.8,E18.8)', i, rOH, energy
+
+        ! Print XYZ geometry
+        write(99, '(I1)') NATOM
+        write(99, *)
+        write(99, '(A,F12.6,F12.6,F12.6)') "O ", x(1,1) / ANG, y(1,1) / ANG, z(1,1) / ANG
+        write(99, '(A,F12.6,F12.6,F12.6)') "H ", x(2,1) / ANG, y(2,1) / ANG, z(2,1) / ANG
+        write(99, '(A,F12.6,F12.6,F12.6)') "H ", x(3,1) / ANG, y(3,1) / ANG, z(3,1) / ANG
+    end do
+
+    ! close file
+    close(99)
+end program

src/h2o_scan_cvrqd.f90

@@ -0,0 +1,67 @@
+! To compile:
+! gfortran -c h2o_scan_cvrqd.f90
+! gfortran -fopenmp -L./ -L../water_potentials/ libabin.a ../water_potentials/libwater.a h2o_scan_cvrqd.o -labin -lwater -lm -lstdc++ -o cvrqd_scan
+program cvrqd_scan
+    use mod_force_h2o
+    use mod_const, only: DP, ANG
+    use mod_init, only: read_xyz_file
+    implicit none
+
+    integer, parameter :: NATOM = 3
+    integer, parameter :: NBEADS = 1
+    integer, parameter :: WATPOT = 1
+
+    ! Displacement for H atom at each timestep
+    real(DP) :: delta = 0.0005_DP
+    integer, parameter :: NSTEPS = 5000
+
+    real(DP) :: x(3, 1), y(3, 1), z(3, 1)
+    real(DP) :: fx(3, 1), fy(3, 1), fz(3, 1)
+    real(DP) :: energy, rOH
+    integer :: i, funit
+    character(len=2) :: atom_names(3)
+    character(len=256) :: fname
+    !character(len=20) :: h2opot = 'cvrqd'
+
+    atom_names = (/ "O", "H", "H" /)
+
+    fname = "water.xyz"
+    open (newunit=funit, file=fname, action="read")
+    call read_xyz_file(funit, fname, atom_names, natom, nbeads, x, y, z)
+    close (funit)
+
+    ! Convert from Angstromgs to Bohr (atomic units)
+    x(:, 1) = x(:, 1) * ANG
+    y(:, 1) = y(:, 1) * ANG
+    z(:, 1) = z(:, 1) * ANG
+
+    ! print for output
+    print*, "Step   ", "rOH   ", "Energy"
+
+    open(99, file="water_geometries.xyz", action="write")
+    ! steps for stretching OH bond
+    do i = 1, NSTEPS
+        ! Calculate the initial OH bond length (rOH) from initial bond positions
+        rOH = sqrt((x(2, 1) - x(1, 1))**2 + (y(2, 1) - y(1, 1))**2 + (z(2, 1) - z(1, 1))**2)
+
+        ! Modify the OH bond length by changing the position of a single hydrogen atom
+        x(2, 1) = x(2, 1) + delta
+
+        ! Call the force_water function to calculate forces and energy
+        energy = 0.0D0
+        call force_h2o(x, y, z, fx, fy, fz, energy, natom, nbeads)
+
+        ! Print the current OH bond length and energy
+        print '(I4,F16.8,E18.8)', i, rOH, energy
+
+        ! Print XYZ geometry
+        write(99, '(I1)') NATOM
+        write(99, *)
+        write(99, '(A,F12.6,F12.6,F12.6)') "O ", x(1,1) / ANG, y(1,1) / ANG, z(1,1) / ANG
+        write(99, '(A,F12.6,F12.6,F12.6)') "H ", x(2,1) / ANG, y(2,1) / ANG, z(2,1) / ANG
+        write(99, '(A,F12.6,F12.6,F12.6)') "H ", x(3,1) / ANG, y(3,1) / ANG, z(3,1) / ANG
+    end do
+
+    ! close file
+    close(99)
+end program

src/h2o_scan_schwenke.f90

@@ -0,0 +1,66 @@
+! To compile:
+! gfortran -c h2o_scan_schwenke.f90
+! gfortran -fopenmp -L./ -L../water_potentials/ libabin.a ../water_potentials/libwater.a h2o_scan_schwenke.o -labin -lwater -lm -lstdc++ -o schwenke_scan
+program schwenke_scan
+    use mod_force_h2o
+    use mod_const, only: DP, ANG
+    use mod_init, only: read_xyz_file
+    implicit none
+
+    integer, parameter :: NATOM = 3
+    integer, parameter :: NBEADS = 1
+    integer, parameter :: WATPOT = 1
+
+    ! Displacement for H atom at each timestep
+    real(DP) :: delta = 0.0005_DP
+    integer, parameter :: NSTEPS = 5000
+
+    real(DP) :: x(3, 1), y(3, 1), z(3, 1)
+    real(DP) :: fx(3, 1), fy(3, 1), fz(3, 1)
+    real(DP) :: energy, rOH
+    integer :: i, funit
+    character(len=2) :: atom_names(3)
+    character(len=256) :: fname
+
+    atom_names = (/ "O", "H", "H" /)
+
+    fname = "water.xyz"
+    open (newunit=funit, file=fname, action="read")
+    call read_xyz_file(funit, fname, atom_names, natom, nbeads, x, y, z)
+    close (funit)
+
+    ! Convert from Angstromgs to Bohr (atomic units)
+    x(:, 1) = x(:, 1) * ANG
+    y(:, 1) = y(:, 1) * ANG
+    z(:, 1) = z(:, 1) * ANG
+
+    ! print for output
+    print*, "Step   ", "rOH   ", "Energy"
+
+    open(99, file="water_geometries.xyz", action="write")
+    ! steps for stretching OH bond
+    do i = 1, NSTEPS
+        ! Calculate the initial OH bond length (rOH) from initial bond positions
+        rOH = sqrt((x(2, 1) - x(1, 1))**2 + (y(2, 1) - y(1, 1))**2 + (z(2, 1) - z(1, 1))**2)
+
+        ! Modify the OH bond length by changing the position of a single hydrogen atom
+        x(2, 1) = x(2, 1) + delta
+
+        ! Call the force_water function to calculate forces and energy
+        energy = 0.0D0
+        call force_h2o(x, y, z, fx, fy, fz, energy, natom, nbeads)
+
+        ! Print the current OH bond length and energy
+        print '(I4,F16.8,E18.8)', i, rOH, energy
+
+        ! Print XYZ geometry
+        write(99, '(I1)') NATOM
+        write(99, *)
+        write(99, '(A,F12.6,F12.6,F12.6)') "O ", x(1,1) / ANG, y(1,1) / ANG, z(1,1) / ANG
+        write(99, '(A,F12.6,F12.6,F12.6)') "H ", x(2,1) / ANG, y(2,1) / ANG, z(2,1) / ANG
+        write(99, '(A,F12.6,F12.6,F12.6)') "H ", x(3,1) / ANG, y(3,1) / ANG, z(3,1) / ANG
+    end do
+
+    ! close file
+    close(99)
+end program