Reference files added for Schwenke testing

src/force_h2o.F90

@@ -70,7 +70,7 @@ subroutine force_h2o_schwenke(x, y, z, fx, fy, fz, Eclas, natom, nbeads)
       integer, intent(in) :: natom, nbeads
       ! Internal water coordinates
       real(DP) :: rOH1, rOH2, aHOH_rad
-      real(DP) :: rij(nbeads, 3)
+      real(DP) :: rij(1, 3)
       real(DP) :: Epot(nbeads)
       integer :: iw
 
@@ -116,26 +116,26 @@ subroutine numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads)
       real(DP) :: z_new_forward(natom, nbeads)
 
       ! This is the energy for the currrent geometry that has already been calculated
-      real(DP), intent(in) :: Epot(nbeads) !nbeads
+      real(DP), intent(in) :: Epot(nbeads)
 
       ! Internal water coordinates
       real(DP) :: new_rOH1, new_rOH2, new_aHOH_rad
       real(DP) :: new_rij(nbeads, natom)
 
       ! Schwenke calculated peterbed geometry energy
-      real(DP) :: Epot_delta(nbeads) !nbeads
+      real(DP) :: Epot_delta(1)
 
       real(DP) :: Eclas_orig, Eclas_plus
       real(DP) :: delta = 5.0E-5_DP
       integer :: i, j, k
 
       ! Calculate forces numerically using central differences
-      do j = 1, nbeads !nbeads
+      do j = 1, nbeads
 
          ! Save the original energy
          Eclas_orig = Epot(j)
 
-         do i = 1, natom !natom
+         do i = 1, natom
 
             do k = 1, 3 ! x, y, z
 
@@ -157,11 +157,11 @@ subroutine numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads)
                ! 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]
+               new_rij(1, :) = [new_rOH1, new_rOH2, new_aHOH_rad]
 
-               call h2o_pot_schwenke(new_rij, Epot_delta, nbeads)
+               call h2o_pot_schwenke(new_rij, Epot_delta(1), nbeads)
 
-               Eclas_plus = Epot_delta(j)
+               Eclas_plus = Epot_delta(1)
 
                ! Calculate the numerical force
                select case (k)

tests/H2O_SCHWENKE/energies.dat.ref

@@ -0,0 +1,2 @@
+ #        Time[fs] E-potential           E-kinetic     E-Total    E-Total-Avg
+           0.48    0.1056085493E-03    0.1328066564E-02    0.1433675113E-02    0.1433675113E-02

tests/H2O_SCHWENKE/forces.xyz.ref

@@ -0,0 +1,5 @@
+           3
+net force: -0.36850E-04 -0.22684E-04  0.21286E-06 torque force:  0.34526E-07  0.32664E-06  0.13444E-04
+O  -0.1084501047E-01  0.4654549435E-02 -0.4100915718E-05
+H   0.8347905491E-02 -0.6517930688E-02  0.5991183457E-05
+H   0.2460254570E-02  0.1840696957E-02 -0.1677405640E-05

tests/H2O_SCHWENKE/input.in

@@ -10,6 +10,7 @@ h2opot='schwenke'
 mdtype='MD',		! classical MD
 dt=20.,			! number of steps and timestep
 nstep=1
+nwritef=1,		! write forces
 /
 
 &nhcopt

tests/H2O_SCHWENKE/movie.xyz.ref

@@ -0,0 +1,5 @@
+           3
+Time step:                   1 Sim. Time [au]          20.00
+O     0.74882117E-03    0.11721115E+00    0.52819796E-03
+H     0.75097634E+00   -0.46201846E+00    0.10479453E-02
+H    -0.75994420E+00   -0.46946219E+00    0.10479453E-02