Remove stubs and comments for Ehrenfest method (#87)

This was never implemented so removing for now.

Author: danielhollas

sample_inputs/input.in.cmd.nab

@@ -7,10 +7,10 @@ Short descriptions of key source files.
 | Path     | Description |
 |----------|-------------|
 | abin.F90        | Main routine. |
-| modules.F90     | Modules containing all basic variables. |
+| modules.F90     | Modules containing all basic global settings. |
 | init.F90        | Big ugly init routine. Reads input parameters, restart files and checks for errors in input. |
 | mdstep.f90      | Routines for propagation of equations of motion (velocity Verlet and RESPA). |
-| forces.f90      | Main routine for getting forces (force_clas) and quantum forces in PIMD (force_quant). |
+| forces.f90      | Main wrapper routine for getting forces (force_clas) and quantum forces in PIMD (force_quant). |
 | force_abin.f90  | Routine that calls ab initio BASH interfaces. |
 | nosehoover.f90  | Nosé-Hoover thermostat. |
 | surfacehop.f90  | Surface Hopping dynamics |
@@ -20,9 +20,8 @@ Short descriptions of key source files.
 | shake.f90       | Constraints using SHAKE algorithm. |
 | analysis.f90    | Driver routine for printing output. |
 | arrays.f90      | Module containing all dynamically allocated arrays.| 
-| random.f90      | Random number generator. |
+| random.f90      | Pseudo-random number generator. |
 | vinit.f90       | Initial velocities. |
-| potentials.f90  | Analytical potentials and hessians for harmonic and Morse oscillators and others. |
-| density.f90     | Evaluates distance, angle and dihedrals distributions during dynamics. | 
+| potentials.f90  | Analytical potentials and hessians (harmonic, Morse etc.)
 | estimators.f90  | Energy and heat capacity estimators for PIMD. | 
 | transform.F90   | Coordinate transformations for PIMD. |

sample_inputs/input.pdb

@@ -62,7 +62,7 @@ program abin
    end if
 
    ! Surface hopping initialization
-   if (ipimd == 2 .or. ipimd == 4) then
+   if (ipimd == 2) then
       call sh_init(x, y, z, vx, vy, vz)
    else if (ipimd == 5 .and. pot == '_tera_') then
       call sh_init(x, y, z, vx, vy, vz)
@@ -117,12 +117,16 @@ program abin
          call force_clas(fxc_diff, fyc_diff, fzc_diff, x, y, z, eclas, pot)
       end if
 
-      ! setting initial values for SURFACE HOPPING
-      if (ipimd == 2 .or. ipimd == 4) then
+      ! setting initial values for surface hopping
+      if (ipimd == 2) then
          do itrj = 1, ntraj
-            if (irest /= 1) call get_nacm(itrj)
+            if (irest /= 1) then
+               call get_nacm(itrj)
+            end if
             call move_vars(vx, vy, vz, vx_old, vy_old, vz_old, itrj)
-            if (pot == '_tera_' .or. restrain_pot == '_tera_') call move_new2old_terash()
+            if (pot == '_tera_' .or. restrain_pot == '_tera_') then
+               call move_new2old_terash()
+            end if
          end do
       else if (ipimd == 5 .and. pot == '_tera_') then
          call move_new2old_terash()

sample_inputs/input.top

@@ -29,10 +29,6 @@ module mod_arrays
    ! for multiple time-stepping and ring contraction ala O. Marsalek
    real(DP), allocatable :: fxc_diff(:, :), fyc_diff(:, :), fzc_diff(:, :)
    real(DP) :: epot_diff = 0.0D0 ! not really array, but for now let's keep it here
-   !ehrenfest forces
-   real(DP), allocatable :: fxeh(:, :), fyeh(:, :), fzeh(:, :) ! interpolation forces
-   real(DP), allocatable :: fxeh_old(:, :), fyeh_old(:, :), fzeh_old(:, :) ! interpolation forces
-   real(DP), allocatable :: fxmid(:, :), fymid(:, :), fzmid(:, :) ! mid step approximative forces
    save
 contains
 
@@ -83,9 +79,6 @@ subroutine allocate_arrays(natomalloc, nwalkalloc)
       vx_old = vx; vy_old = vx_old; vz_old = vx_old
       transxv = vx; transyv = transxv; transzv = transxv
 
-      ! TODO-EH: we need to allocate fxc according to the number of states
-      ! but this subroutine is called before we now the number of states
-      ! so leave this as is and reallocate elsewhere
       allocate (fxc(natomalloc, nwalkalloc))
       allocate (fyc(natomalloc, nwalkalloc))
       allocate (fzc(natomalloc, nwalkalloc))
@@ -114,30 +107,6 @@ subroutine allocate_arrays(natomalloc, nwalkalloc)
 
    end subroutine allocate_arrays
 
-   !EHRENFEST gradients and midstep approximative forces for all states
-   subroutine allocate_ehrenfest(natom, nstate)
-      integer, intent(in) :: nstate
-      integer, intent(in) :: natom
-      allocate (fxeh(natom, nstate))
-      allocate (fyeh(natom, nstate))
-      allocate (fzeh(natom, nstate))
-      allocate (fxeh_old(natom, nstate))
-      allocate (fyeh_old(natom, nstate))
-      allocate (fzeh_old(natom, nstate))
-      allocate (fxmid(natom, nstate))
-      allocate (fymid(natom, nstate))
-      allocate (fzmid(natom, nstate))
-      fxeh = 0.0D0
-      fyeh = 0.0D0
-      fzeh = 0.0D0
-      fxeh_old = 0.0D0
-      fyeh_old = 0.0D0
-      fzeh_old = 0.0D0
-      fxmid = 0.0D0
-      fymid = 0.0D0
-      fzmid = 0.0D0
-   end subroutine allocate_ehrenfest
-
    subroutine deallocate_arrays
       ! If x is allocated, all of them should be.
       ! If not, something is horribly wrong anyway
@@ -159,20 +128,6 @@ subroutine deallocate_arrays
       if (allocated(fxc_diff)) then
          deallocate (fxc_diff); deallocate (fyc_diff); deallocate (fzc_diff); 
       end if
-
-      ! Ehrenfest variables (hmm, are these really needed??
-      ! Would be better to just use the existing ones...
-      if (allocated(fxeh)) then
-         deallocate (fxeh)
-         deallocate (fyeh)
-         deallocate (fzeh)
-         deallocate (fxeh_old)
-         deallocate (fyeh_old)
-         deallocate (fzeh_old)
-         deallocate (fxmid)
-         deallocate (fymid)
-         deallocate (fzmid)
-      end if
    end subroutine deallocate_arrays
 
 end module mod_arrays

src/README.md

@@ -61,15 +61,14 @@ subroutine force_abin(x, y, z, fx, fy, fz, eclas, chpot, walkmax)
       end do
       close (unit=MAXUNITS + iw)
 
-!     Surface hopping or Ehrenfest
-      if (ipimd == 2 .or. ipimd == 4) then
+      if (ipimd == 2) then
 
          open (unit=MAXUNITS + iw + 2 * walkmax, file='state.dat')
          write (MAXUNITS + iw + 2 * walkmax, '(I2)') nstate
 
-!        Diagonal of tocalc holds info about needed forces
-!        tocalc(x,x)= 1 -> compute forces for electronic state X
-!        totalc for Ehrenfest muset be set just for required states according to c coef. TO-DO in ehrenfest enrehfest_forces
+         ! Diagonal of tocalc holds info about needed forces
+         ! tocalc(x,x) = 1 -> compute forces for electronic state X
+         ! off-diagonal elements correspond to non-adiabatic couplings
          do ist1 = 1, nstate
             write (MAXUNITS + iw + 2 * walkmax, '(I1,A1)', advance='no') tocalc(ist1, ist1), ' '
          end do
@@ -82,10 +81,7 @@ subroutine force_abin(x, y, z, fx, fy, fz, eclas, chpot, walkmax)
          open (unit=MAXUNITS + iw + 2 * walkmax, file='state.dat')
          write (MAXUNITS + iw + 2 * walkmax, '(I2)') nstate_lz !How many el. states
 
-         !do ist1=1,nstate_lz
-         !   write(MAXUNITS+iw+2*walkmax,'(I1,A1)',advance='no')tocalc_lz(ist1),' '
-         !end do
-         !First we have singlets, then triplets
+         ! First we have singlets, then triplets
          do ist1 = 1, nstate_lz
             if (tocalc_lz(ist1) == 1) write (MAXUNITS + iw + 2 * walkmax, '(I2,A1)') ist1, ' ' !Number of gradient state
          end do
@@ -114,7 +110,7 @@ subroutine force_abin(x, y, z, fx, fy, fz, eclas, chpot, walkmax)
 
 !     for SH, pass the 4th parameter: precision of forces as 10^(-force_accu1)
 !     TODO: This should not be hard-coded
-      if (ipimd == 2 .or. ipimd == 4 .or. ipimd == 5) then
+      if (ipimd == 2 .or. ipimd == 5) then
          write (chsystem, '(A60,I3,A12)') chsystem, 7, ' < state.dat'
       end if
 
@@ -163,12 +159,11 @@ subroutine force_abin(x, y, z, fx, fy, fz, eclas, chpot, walkmax)
       eclas = eclas + temp1
 ! SH
       ! TODO: Have each state in different file?
-      if (ipimd == 2 .or. ipimd == 4) then
+      if (ipimd == 2) then
          en_array(1, iw) = temp1
          do ist1 = 2, nstate
             read (MAXUNITS + iw, *) en_array(ist1, iw)
          end do
-         ! TODO-EH: eclas must be correctly overwritten later
          eclas = en_array(istate(iw), iw)
       end if
 ! LZ
@@ -189,11 +184,7 @@ subroutine force_abin(x, y, z, fx, fy, fz, eclas, chpot, walkmax)
          eclas = en_array_lz(istate_lz, 1)
       end if
 
-!     TODO-EH: Read additional forces probably somewhere here, use second index (iw) for different states
-!     Actually, we should make a routine read_engrad() and make it general for both EH and SH
-!     always read energies, read forces based on tocalc
-
-      if (ipimd == 2 .or. ipimd == 4) then
+      if (ipimd == 2) then
          iost = read_forces(fx, fy, fz, natqm, tocalc(istate(iw), istate(iw)), MAXUNITS + iw)
       else if (ipimd == 5) then
          iost = read_forces(fx, fy, fz, natqm, tocalc_lz(istate_lz), MAXUNITS + iw) !Save only the computed state

src/abin.F90

@@ -198,7 +198,7 @@ subroutine force_wrapper(x, y, z, fx, fy, fz, e_pot, chpot, walkmax)
    case ("_cp2k_")
       call force_cp2k(x, y, z, fx, fy, fz, eclas, walkmax)
    case ("_tera_")
-      if (ipimd == 2 .or. ipimd == 4 .or. ipimd == 5) then
+      if (ipimd == 2 .or. ipimd == 5) then
          call force_terash(x, y, z, fx, fy, fz, eclas)
       else
          call force_tera(x, y, z, fx, fy, fz, eclas, walkmax)

src/arrays.F90

@@ -16,6 +16,7 @@ subroutine init(dt)
    use mod_const
    use mod_interfaces, only: print_compile_info, omp_set_num_threads, print_runtime_info
    use mod_cmdline, only: get_cmdline
+   use mod_error, only: fatal_error
    use mod_files
    use mod_arrays
    use mod_array_size
@@ -222,10 +223,10 @@ subroutine init(dt)
          ipimd = 2
       case ('minimization')
          ipimd = 3
-      case ('ehrenfest')
-         ipimd = 4
       case ('landau_zener')
          ipimd = 5
+      case default
+         call fatal_error(__FILE__, __LINE__, 'invalid mdtype in '//trim(chinput))
       end select
    end if
 
@@ -265,8 +266,6 @@ subroutine init(dt)
          print '(a)', '           Surface Hopping MD                 '
       case (3)
          print '(a)', '              Minimization                    '
-      case (4)
-         print '(a)', '              Ehrenfest MD                    '
       case (5)
          print '(a)', '             Landau Zener MD                  '
       end select
@@ -348,9 +347,6 @@ subroutine init(dt)
 
 !  allocate all basic arrays and set them to 0.0d0
    call allocate_arrays(natom, nwalk + 1)
-!  Ehrenfest require larger array since gradients for all of the states are need
-!  TODO: We should really make this differently..
-   if (ipimd == 4) call allocate_ehrenfest(natom, nstate)
 
    if (iplumed == 1) then
       call plumed_init()
@@ -417,7 +413,7 @@ subroutine init(dt)
    read (150, nhcopt)
    rewind (150)
 
-   if (ipimd == 2 .or. ipimd == 4) then
+   if (ipimd == 2) then
       read (150, sh)
       rewind (150)
       integ = tolower(integ)
@@ -445,15 +441,15 @@ subroutine init(dt)
 
    if (pot == '_tera_' .or. restrain_pot == '_tera_') then
       call initialize_tc_servers()
-      if (ipimd == 2 .or. ipimd == 4 .or. ipimd == 5) then
+      if (ipimd == 2 .or. ipimd == 5) then
          call init_terash(x, y, z)
       end if
    end if
 
-!-----HERE WE CHECK FOR ERRORS IN INPUT-----------
+   ! Check whether input parameters make sense
    call check_inputsanity()
 
-!     resetting number of walkers to 1 in case of classical simulation
+   ! resetting number of walkers to 1 in case of classical simulation
    if (ipimd == 0) then
       if (my_rank == 0) then
          write (*, *) 'Using velocity Verlet integrator'
@@ -465,8 +461,9 @@ subroutine init(dt)
       ! algorithm as well
    end if
 
-!for surface hopping and ehrenfest
-   if (ipimd == 2 .or. ipimd == 4 .or. ipimd == 5) then
+   ! Selecting proper integrator for a given MD type
+   ! (controlled by variable 'md')
+   if (ipimd == 2 .or. ipimd == 5) then
       nwalk = ntraj !currently 1
       md = 2 ! velocity verlet
       nabin = 1
@@ -648,7 +645,7 @@ subroutine init(dt)
          write (*, nml=nhcopt, delim='APOSTROPHE')
          write (*, *)
       end if
-      if (ipimd == 2 .or. ipimd == 4) then
+      if (ipimd == 2) then
          write (*, nml=sh, delim='APOSTROPHE')
          write (*, *)
       end if
@@ -911,7 +908,7 @@ subroutine check_inputsanity()
          write (*, *) chknow
          if (iknow /= 1) error = 1
       end if
-      if (inose == 1 .and. (ipimd == 2 .or. ipimd == 4) .and. en_restraint == 0) then
+      if (inose == 1 .and. (ipimd == 2) .and. en_restraint == 0) then
          write (*, *) 'Thermostating is not meaningful for surface hopping simulation.'
          write (*, *) chknow
          if (iknow /= 1) error = 1

src/force_abin.F90

@@ -80,11 +80,9 @@ end subroutine thermostat
 !  Contains propagation of normal modes according to:
 !  eq 23 from J. Chem. Phys. 133, 124104 2010
    subroutine verletstep(x, y, z, px, py, pz, amt, dt, eclas, fxc, fyc, fzc)
-      use mod_general, only: pot, ipimd, inormalmodes, en_restraint
+      use mod_general, only: pot, inormalmodes, en_restraint
       use mod_nhc, only: inose
       use mod_interfaces, only: force_clas, propagate_nm
-      ! Not yet implemented
-      !use mod_sh, only: ehrenfest_forces
       use mod_en_restraint
       real(DP), intent(inout) :: x(:, :), y(:, :), z(:, :)
       real(DP), intent(inout) :: fxc(:, :), fyc(:, :), fzc(:, :)
@@ -107,13 +105,6 @@ subroutine verletstep(x, y, z, px, py, pz, amt, dt, eclas, fxc, fyc, fzc)
 
       call force_clas(fxc, fyc, fzc, x, y, z, eclas, pot)
 
-      if (ipimd == 4) then
-         ! This is only a stub for now
-         write (*, *) 'ERROR: Ehrenfest MD not implemented yet!!'
-         call abinerror('verletstep')
-         !call ehrenfest_forces(x, y, z, fxc, fyc, fzc, px, py, pz, dt, eclas)
-      end if
-
       call shiftP(px, py, pz, fxc, fyc, fzc, dt / 2)
 
       if (inose > 0) call thermostat(px, py, pz, amt, dt / 2)