NaI analytic potential for FSSH (#178)

2-state analytic potential for NaI implemented for surface hopping dynamics. The original potential is diabatic but was adiabatized in the code. 3ps of dynamics can be done in 8 seconds when one is not writing too much on the disk. Accessing the disk is the most demanding part of the simulation. The implementation was compared to the exact QD.
PHOTOX · Jul 3, 2024 · 46e61e5 · 46e61e5
1 parent 675ae72
commit 46e61e5
Show file tree

Hide file tree

Showing 16 changed files with 1,574 additions and 6 deletions.
diff --git a/src/Makefile b/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 h2o_schwenke.o h2o_cvrqd.o force_h2o.o force_cp2k.o sh_integ.o surfacehop.o landau_zener.o\
+          remd.o force_bound.o water.o h2o_schwenke.o h2o_cvrqd.o force_h2o.o force_cp2k.o sh_integ.o surfacehop.o landau_zener.o potentials_sh.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
 

diff --git a/src/forces.F90 b/src/forces.F90
@@ -133,6 +133,7 @@ subroutine force_wrapper(x, y, z, fx, fy, fz, e_pot, chpot, walkmax)
    use mod_sbc, only: force_sbc, isbc
    use mod_plumed, only: iplumed, force_plumed
    use mod_potentials, only: force_harmonic_rotor, force_harmonic_oscillator, force_morse, force_doublewell
+   use mod_potentials_sh, only: force_nai
    use mod_splined_grid, only: force_splined_grid
    use mod_cp2k, only: force_cp2k
    use mod_shell_interface, only: force_abin
@@ -173,6 +174,8 @@ subroutine force_wrapper(x, y, z, fx, fy, fz, e_pot, chpot, walkmax)
       call force_cp2k(x, y, z, fx, fy, fz, eclas, walkmax)
    case ("_tcpb_")
       call force_tcpb(x, y, z, fx, fy, fz, eclas, walkmax)
+   case ("_nai_")
+      call force_nai(x, y, z, fx, fy, fz, eclas)
    case ("_tera_")
       if (ipimd == 2 .or. ipimd == 5) then
          call force_terash(x, y, z, fx, fy, fz, eclas)

diff --git a/src/init.F90 b/src/init.F90
@@ -40,7 +40,8 @@ subroutine init(dt)
       use mod_nhc
       use mod_estimators
       use mod_potentials
-      use mod_sh_integ, only: phase
+      use mod_potentials_sh
+      use mod_sh_integ, only: phase, nstate
       use mod_sh, only: read_sh_input, print_sh_input
       use mod_lz
       use mod_qmmm, only: natqm, natmm
@@ -518,6 +519,9 @@ subroutine init(dt)
       if (pot == '_mm_' .or. pot_ref == '_mm_') then
          call initialize_mm(natom, atnames, mm_types, q, LJ_rmin, LJ_eps)
       end if
+      if (pot == '_nai_' .or. pot_ref == '_nai_') then
+         call nai_init(natom, nwalk, ipimd, nstate)
+      end if
 
       if (my_rank == 0) then
          call print_simulation_parameters()

diff --git a/src/potentials_sh.F90 b/src/potentials_sh.F90
@@ -0,0 +1,170 @@
+! Various analytical potentials for surface hopping:
+! - NaI potential and couplings in adiabatic basis
+
+! Created by Jiri Janos
+
+module mod_potentials_sh
+   use mod_const, only: DP
+   implicit none
+   public
+   private :: nai
+
+   ! We use derived types to encapsulate potential parameters
+   ! https://fortran-lang.org/learn/quickstart/derived_types
+
+   ! Parameters for the NaI model as defined in the following articles
+   ! https://doi.org/10.1063/1.4919780
+   ! https://doi.org/10.1063/1.456377
+   ! These articles contain diabatic basis while here we use adiabatic. In this code, we first calculate the diabatic quantities and
+   ! then we convert them to adiabatic using a simple diagonalization of 2 states. The equations for diabatic quantites can be found
+   ! in the code or in appendix of this article: https://doi.org/10.1063/5.0071376 (sign is not correct) or equation 1 of this paper
+   ! https://doi.org/10.1063/1.1377030.
+   type :: nai_params
+      ! Diabatic state VX parameters
+      real(DP) :: a2 = 2760D0
+      real(DP) :: b2 = 2.398D0
+      real(DP) :: c2 = 11.3D0
+      real(DP) :: lp = 0.408D0
+      real(DP) :: lm = 6.431D0
+      real(DP) :: rho = 0.3489D0
+      real(DP) :: de0 = 2.075D0
+      real(DP) :: e2 = 14.399613877582553D0 ! elementary charge in eV/angs units
+      ! Diabatic state VA parameters
+      real(DP) :: a1 = 0.813D0
+      real(DP) :: beta1 = 4.08D0
+      real(DP) :: r0 = 2.67D0
+      ! Diabatic coupling VXA parameters
+      real(DP) :: a12 = 0.055D0
+      real(DP) :: beta12 = 0.6931D0
+      real(DP) :: rx = 6.93D0
+   end type
+   type(nai_params) :: nai
+   save
+contains
+
+   ! NaI potential
+   subroutine nai_init(natom, nwalk, ipimd, nstate)
+      use mod_error, only: fatal_error
+      use mod_utils, only: real_positive
+      integer, intent(in) :: natom, nwalk, ipimd, nstate
+
+      if (natom /= 2) then
+         call fatal_error(__FILE__, __LINE__, &
+            & 'NaI potential is only for 2 particles.')
+      end if
+
+      if (nwalk /= 1) then
+         call fatal_error(__FILE__, __LINE__, &
+            & 'NaI model requires only 1 walker.')
+      end if
+
+      if (ipimd /= 2) then
+         call fatal_error(__FILE__, __LINE__, &
+            & 'NaI model works only with SH dynamics.')
+      end if
+
+      if (nstate /= 2) then
+         call fatal_error(__FILE__, __LINE__, &
+            & 'NaI model is implemented for 2 states only.')
+      end if
+
+   end subroutine nai_init
+
+   subroutine force_nai(x, y, z, fx, fy, fz, eclas)
+      use mod_sh, only: en_array, istate, nacx, nacy, nacz
+      use mod_const, only: ANG, AUTOEV
+      real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
+      real(DP), intent(out) :: fx(:, :), fy(:, :), fz(:, :)
+      real(DP), intent(out) :: eclas
+      real(DP) :: VX, VA, VXA, dVX, dVA, dVXA ! diabatic hamiltonian and its derivatives
+      real(DP) :: E1, E2, d12, dE1, dE2 ! adiabatic hamiltonian and derivatives of energies
+      real(DP) :: r, fr, dx, dy, dz
+
+      ! NOTE: The original potential is defined in diabatic basis. For ABIN's purposes, we need to transfer to adiabatic basis,
+      ! which can be easily done for two states. However, the adiabatic formulas are very long and complicated to both code and
+      ! read. Therefore, we calculate the diabatic quantities and their derivatives and then construct the adiabatic states and
+      ! couplings from diabatic ones. It should be easier to read and also more efficient. Note that the diabatic potential is
+      ! in eV and Angstrom so conversions are necessary.
+
+      ! distance between atoms
+      dx = x(2, 1) - x(1, 1)
+      dy = y(2, 1) - y(1, 1)
+      dz = z(2, 1) - z(1, 1)
+      r = dx**2 + dy**2 + dz**2
+      r = dsqrt(r)
+
+      ! normalized distance
+      dx = dx / r
+      dy = dy / r
+      dz = dz / r
+
+      ! converting distance to Angstrom as the potential was done
+      r = r / ANG
+
+      ! calculating diabatic hamiltonian
+      VX = (nai%a2 + (nai%b2 / r)**8) * dexp(-r / nai%rho) - nai%e2 / r - nai%e2 * (nai%lp + nai%lm) / 2 / r**4 - &
+           nai%c2 / r**6 - 2 * nai%e2 * nai%lm * nai%lp / r**7 + nai%de0
+      VA = nai%a1 * dexp(-nai%beta1 * (r - nai%r0))
+      VXA = nai%a12 * dexp(-nai%beta12 * (r - nai%rx)**2)
+      ! calculating derivatives of diabatic hamiltonian
+      dVX = -(nai%a2 + (nai%b2 / r)**8) * dexp(-r / nai%rho) / nai%rho - 8.0D0 * nai%b2**8 / r**9 * dexp(-r / nai%rho) + &
+            14.0D0 * nai%e2 * nai%lm * nai%lp / r**8 + 6.0D0 * nai%c2 / r**7 + 2.0D0 * nai%e2 * (nai%lp + nai%lm) / r**5 + &
+            nai%e2 / r**2
+      dVA = -nai%beta1 * VA
+      dVXA = -2.0D0 * nai%beta12 * (r - nai%rx) * VXA
+      ! converting them to a.u. as they are in eV or eV/Angstrom
+      VX = VX / AUTOEV
+      VA = VA / AUTOEV
+      VXA = VXA / AUTOEV
+      dVX = dVX / AUTOEV / ANG
+      dVA = dVA / AUTOEV / ANG
+      dVXA = dVXA / AUTOEV / ANG
+
+      ! adiabatic potentials
+      E1 = (VA + VX) / 2.0D0 - dsqrt((VA - VX)**2.0D0 + 4.0D0 * VXA**2.0D0) / 2.0D0
+      E2 = (VA + VX) / 2.0D0 + dsqrt((VA - VX)**2.0D0 + 4.0D0 * VXA**2.0D0) / 2.0D0
+      ! nonadiabatic coupling vector in the reduced system
+      d12 = -(VXA * (dVA - dVX) + (-VA + VX) * dVXA) / (VA**2.0D0 - 2.0D0 * VA * VX + VX**2.0D0 + 4.0D0 * VXA**2.0D0)
+      d12 = d12 / ANG ! one more conversion necessary
+      ! derivatives of energies
+      dE1 = (dVA + dVX) / 2.0D0 - (2.0D0 * (VA - VX) * (dVA - dVX) + 8.0D0 * VXA * dVXA) / &
+            (4.0D0 * dsqrt((VA - VX)**2.0D0 + 4.0D0 * VXA**2.0D0))
+      dE2 = (dVA + dVX) / 2.0D0 + (2.0D0 * (VA - VX) * (dVA - dVX) + 8.0D0 * VXA * dVXA) / &
+            (4.0D0 * dsqrt((VA - VX)**2.0D0 + 4.0D0 * VXA**2.0D0))
+
+      ! saving electronic energies for SH
+      en_array(1) = E1
+      en_array(2) = E2
+
+      ! saving classical energy for Verlet
+      eclas = en_array(istate)
+
+      ! calculate force (-dE/dr) on the propagated state
+      if (istate == 1) fr = -dE1
+      if (istate == 2) fr = -dE2
+
+      ! projecting the force on the cartesian coordinates
+      fx(1, 1) = -fr * dx
+      fx(2, 1) = -fx(1, 1)
+      fy(1, 1) = -fr * dy
+      fy(2, 1) = -fy(1, 1)
+      fz(1, 1) = -fr * dz
+      fz(2, 1) = -fz(1, 1)
+
+      ! projecting coupling to the cartesian coordinates
+      nacx(1, 1, 2) = -d12 * dx
+      nacy(1, 1, 2) = -d12 * dy
+      nacz(1, 1, 2) = -d12 * dz
+      nacx(1, 2, 1) = -nacx(1, 1, 2)
+      nacy(1, 2, 1) = -nacy(1, 1, 2)
+      nacz(1, 2, 1) = -nacz(1, 1, 2)
+      nacx(2, 1, 2) = d12 * dx
+      nacy(2, 1, 2) = d12 * dy
+      nacz(2, 1, 2) = d12 * dz
+      nacx(2, 2, 1) = -nacx(2, 1, 2)
+      nacy(2, 2, 1) = -nacy(2, 1, 2)
+      nacz(2, 2, 1) = -nacz(2, 1, 2)
+
+   end subroutine force_nai
+
+end module mod_potentials_sh
diff --git a/src/surfacehop.F90 b/src/surfacehop.F90
@@ -287,6 +287,7 @@ subroutine get_nacm(pot)
 
       ! In TeraChem SH interface, we already got NACME
       if (pot == '_tera_') return
+      if (pot == '_nai_') return ! for NaI model, the couplings are calcualted together with forces
 
       ! Check whether we need to compute any NACME
       num_nacme = 0
@@ -611,16 +612,14 @@ subroutine surfacehop(x, y, z, vx, vy, vz, vx_old, vy_old, vz_old, dt, eclas)
 
       ! First, calculate NACME
       if (inac == 0) then
-
          ! For TeraChem MPI / FMS interface, NAC are already computed!
-         if (pot /= '_tera_') then
+         if (pot /= '_tera_' .and. pot /= '_nai_') then
             nacx = 0.0D0
             nacy = 0.0D0
             nacz = 0.0D0
             ! Compute and read NACME (MOLPRO-SH interface)
             call get_nacm(pot)
          end if
-
          ! TODO: Should we call this with TeraChem?
          ! I think TC already phases the couplings internally.
          call phase_nacme(nacx_old, nacy_old, nacz_old, nacx, nacy, nacz)

diff --git a/tests/SH_NaI/PES.dat.ref b/tests/SH_NaI/PES.dat.ref
@@ -0,0 +1,6 @@
+ #    Time[fs] Potential energies
+           0.06   -0.8731544800E-03    0.3904840511E-02
+           0.12   -0.8738840794E-03    0.3902676524E-02
+           0.18   -0.8746152732E-03    0.3900509975E-02
+           0.24   -0.8753480644E-03    0.3898340868E-02
+           0.30   -0.8760824558E-03    0.3896169203E-02
diff --git a/tests/SH_NaI/distances.dat.ref b/tests/SH_NaI/distances.dat.ref
@@ -0,0 +1,6 @@
+ # Distances [Angstrom]
+           0.06   7.2911701
+           0.12   7.2908894
+           0.18   7.2906084
+           0.24   7.2903270
+           0.30   7.2900452
diff --git a/tests/SH_NaI/dotprod.dat.ref b/tests/SH_NaI/dotprod.dat.ref
@@ -0,0 +1,6 @@
+ #    Time[fs] dotproduct(i,j) [i=1,nstate-1;j=i+1,nstate]
+           0.06   -0.5679102843E-04 
+           0.12   -0.5690157880E-04 
+           0.18   -0.5701226238E-04 
+           0.24   -0.5712307933E-04 
+           0.30   -0.5723402981E-04 
diff --git a/tests/SH_NaI/input.in b/tests/SH_NaI/input.in
@@ -0,0 +1,37 @@
+This is a sample input file for Surface-Hopping simulation in ABIN 
+
+&general
+pot='_nai_'		   ! where do we obtain forces?
+irest=0,		   ! should we restart from restart.xyz?
+irandom=347110445, ! random seed
+
+mdtype='sh',	! sh = Surface Hopping non adiabatic MD
+nstep=5,    	! number of steps, originally in QD 260000
+dt=2.5,			! timestep [au], originally in QD 0.25
+
+nwrite=1,		! how often some output should be printed (energies, surface hopping probabilities etc.)
+nwritex=1,		! how often to print coordinates?
+nwritev=1,		! how often to print coordinates?
+narchive=90000,	! how often to print coordinates?
+nrest=200,		! how often to print restart files?
+/
+
+&nhcopt
+inose=0,		! NVE ensemble (thermostat turned OFF)
+!temp=0.00,		! Usually, you would read initial velocities from previous ground state sampling (-v option)
+/
+
+&sh
+istate_init=2,		! initial electronic state (1 is ground state)
+nstate=2,		! number of electronic states
+deltaE=2.0,		! maximum energy difference (eV) between states for which we calculate NA coupling
+PopThr=0.00001,         ! minimum population of either state, for which we compute NA coupling
+EnergyDifThr=0.01,      ! maximum energy difference between two consecutive steps
+EnergyDriftThr=0.01,    ! maximum energy drift from initial total energy
+/
+
+&system
+ndist=1,
+dist1=1,
+dist2=2,
+/
diff --git a/tests/SH_NaI/mini.xyz b/tests/SH_NaI/mini.xyz
@@ -0,0 +1,4 @@
+           2
+Time step:                2700 Sim. Time [au]        6750.00
+I    -0.71520124E+00    0.00000000E+00    0.00000000E+00
+Na    0.65762491E+01    0.00000000E+00    0.00000000E+00
diff --git a/tests/SH_NaI/nacm_all.dat.ref b/tests/SH_NaI/nacm_all.dat.ref
@@ -0,0 +1,20 @@
+ Time step:           1
+ NACME between states:           1           2
+   -0.2678921537E+00   -0.0000000000E+00   -0.0000000000E+00
+    0.2678921537E+00    0.0000000000E+00    0.0000000000E+00
+ Time step:           2
+ NACME between states:           1           2
+   -0.2680501091E+00   -0.0000000000E+00   -0.0000000000E+00
+    0.2680501091E+00    0.0000000000E+00    0.0000000000E+00
+ Time step:           3
+ NACME between states:           1           2
+   -0.2682083260E+00   -0.0000000000E+00   -0.0000000000E+00
+    0.2682083260E+00    0.0000000000E+00    0.0000000000E+00
+ Time step:           4
+ NACME between states:           1           2
+   -0.2683668044E+00   -0.0000000000E+00   -0.0000000000E+00
+    0.2683668044E+00    0.0000000000E+00    0.0000000000E+00
+ Time step:           5
+ NACME between states:           1           2
+   -0.2685255443E+00   -0.0000000000E+00   -0.0000000000E+00
+    0.2685255443E+00    0.0000000000E+00    0.0000000000E+00
diff --git a/tests/SH_NaI/pop.dat.ref b/tests/SH_NaI/pop.dat.ref
@@ -0,0 +1,6 @@
+ #    Time[fs] CurrentState   Populations Sum-of-Populations
+           0.06  2   0.00000   1.00000 0.9999999
+           0.12  2   0.00000   1.00000 0.9999999
+           0.18  2   0.00000   1.00000 0.9999999
+           0.24  2   0.00000   1.00000 0.9999999
+           0.30  2   0.00000   1.00000 0.9999999
diff --git a/tests/SH_NaI/prob.dat.ref b/tests/SH_NaI/prob.dat.ref
@@ -0,0 +1,6 @@
+ #    Time[fs] CurrentState   Probabilities
+           0.06  2   0.00000   0.00000
+           0.12  2   0.00000   0.00000
+           0.18  2   0.00000   0.00000
+           0.24  2   0.00000   0.00000
+           0.30  2   0.00000   0.00000