mirror of
https://github.com/nwchemgit/nwchem.git
synced 2026-07-28 22:25:48 -04:00
same as 27006 from 6.5: several enhancements to existing experimental qmd module (from niri)
This commit is contained in:
Edoardo Apra
parent
b9d2753252
commit
e80a6bb61b
11 changed files with 494 additions and 139 deletions
|
|
@ -4,7 +4,10 @@
|
|||
OBJ_OPTIMIZE = qmd_init.o qmd_driver.o qmd_main.o qmd_input.o\
|
||||
qmd_start.o task_qmd.o qmd_nucl_dt2.o qmd_nucl_dt.o\
|
||||
qmd_gradient.o qmd_geom_store.o qmd_thermostat.o\
|
||||
qmd_kinetic.o qmd_write_trj.o qmd_write_rst.o
|
||||
qmd_kinetic.o qmd_write_trj.o qmd_write_rst.o\
|
||||
qmd_ran1.o qmd_gasdev.o qmd_randomize.o\
|
||||
qmd_nucl_pos.o qmd_nucl_vel.o qmd_read_rst.o\
|
||||
qmd_gamdev.o qmd_sumno.o qmd_fix_momentum.o
|
||||
|
||||
LIBRARY=libqmd.a
|
||||
|
||||
|
|
|
|||
|
|
@ -7,12 +7,8 @@ c
|
|||
common/qmd_int/lgeom,igeom,jdebug
|
||||
c
|
||||
c constants
|
||||
double precision kb,au2fs
|
||||
common/qmd_consts/kb,au2fs
|
||||
c
|
||||
c berendsen thermostat
|
||||
logical do_berendsen
|
||||
double precision tau
|
||||
double precision kb,au2fs,amu2au,au2ang
|
||||
common/qmd_consts/kb,au2fs,amu2au,au2ang
|
||||
c
|
||||
c nose-hoover thermostat parameters
|
||||
integer maxchain
|
||||
|
|
@ -22,9 +18,27 @@ c nose-hoover thermostat parameters
|
|||
double precision v_nh(maxchain) ! chain velocities
|
||||
double precision m_nh(maxchain) ! chain masses
|
||||
double precision g_nh(maxchain) ! bath couplings
|
||||
c
|
||||
c berendsen thermostat
|
||||
logical do_berendsen
|
||||
double precision tau
|
||||
c
|
||||
c langevin thermostat
|
||||
logical do_langevin
|
||||
double precision friction
|
||||
integer idum
|
||||
c
|
||||
c velocity rescaling
|
||||
logical do_rescale
|
||||
c
|
||||
c BDP rescaling
|
||||
logical do_BDP
|
||||
c
|
||||
common/qmd_tstats/do_berendsen,tau,
|
||||
& do_nosehoover,r_nh,v_nh,m_nh,g_nh
|
||||
common/qmd_thermostats/do_berendsen,tau,
|
||||
& do_nosehoover,r_nh,v_nh,m_nh,g_nh,
|
||||
& do_langevin,friction,
|
||||
& do_rescale,idum,
|
||||
& do_BDP
|
||||
c
|
||||
c temperature
|
||||
double precision targ_temp
|
||||
|
|
|
|||
|
|
@ -34,13 +34,18 @@ c
|
|||
integer nat ! number of atoms
|
||||
integer nactive ! number of active atoms
|
||||
integer nstep_nucl ! total nuclear steps
|
||||
integer sstep_nucl ! first nuclear step
|
||||
integer istep_nucl ! nuclear md step
|
||||
integer fstep ! how often to remove trans/rot
|
||||
double precision dt_nucl ! nuclear md time step
|
||||
double precision eltime ! elapsed time (in fs)
|
||||
double precision esys ! system energy (from task_gradient)
|
||||
double precision ekin ! kinetic energy from the md part
|
||||
double precision etotal ! total energy: esys + ekin
|
||||
c
|
||||
logical status
|
||||
logical rst_exists
|
||||
logical do_none
|
||||
c
|
||||
character*32 pname
|
||||
integer xyzfileunit
|
||||
|
|
@ -73,7 +78,7 @@ c get rstfilename
|
|||
call util_file_name('qmdrst',.false.,.false.,rstfilename)
|
||||
c
|
||||
c initialize parameters
|
||||
call qmd_init(rtdb,nstep_nucl,dt_nucl)
|
||||
call qmd_init(rtdb,nstep_nucl,dt_nucl,do_none,fstep)
|
||||
c
|
||||
c create geometry object
|
||||
if (.not. geom_create(geom,'geometry'))
|
||||
|
|
@ -106,7 +111,6 @@ c
|
|||
if(.not.ma_push_get(mt_dbl,3*nat,'grads',l_g,k_g)) ! grads
|
||||
$ call errquit(pname//'ma_push_get grads',3*nat,MA_ERR)
|
||||
call dfill(3*nat, 0.0d0, dbl_mb(k_g), 1)
|
||||
c
|
||||
c get the active atoms (also need a check for linearity)
|
||||
c need to implement constraints
|
||||
if (.not. ma_push_get(mt_log,nat,'active atoms',l_act,k_act)) ! active
|
||||
|
|
@ -123,25 +127,68 @@ c get the masses, initial coordinates
|
|||
$ call errquit(pname//'geom_masses_get',geom,GEOM_ERR)
|
||||
if (.not.geom_cart_coords_get(geom,dbl_mb(k_r)))
|
||||
$ call errquit(pname//'geom_cart_coords_get',0,GEOM_ERR)
|
||||
c
|
||||
inquire(file=rstfilename,exist=rst_exists)
|
||||
if (rst_exists) then
|
||||
if (ga_nodeid().eq.0) then
|
||||
write(luout,*) 'Restart file found'
|
||||
write(luout,*) 'Reading positions, velocities, KE,'
|
||||
write(luout,*) 'and last time step from file.'
|
||||
end if
|
||||
c read velocities from file
|
||||
call qmd_read_rst(rstfilename,rstfileunit,nat,dbl_mb(k_v),
|
||||
$ dbl_mb(k_r),ekin,sstep_nucl)
|
||||
c remove global translations and rotations
|
||||
c kinetic energy is the same on exit as was on entry
|
||||
call qmd_fix_momentum(nat,dbl_mb(k_m),dbl_mb(k_v),dbl_mb(k_r),
|
||||
$ ekin,.true.)
|
||||
c update geom and store in rtdb
|
||||
call qmd_geom_store(rtdb,geom,nat,dbl_mb(k_r))
|
||||
else
|
||||
if (ga_nodeid().eq.0) then
|
||||
write(luout,*) 'No restart file found'
|
||||
write(luout,*) 'Beginning with random velocities'
|
||||
end if
|
||||
c randomize velocities for the given temperature
|
||||
call qmd_randomize(nat,dbl_mb(k_m),dbl_mb(k_v))
|
||||
sstep_nucl=1
|
||||
c remove global translations and rotations
|
||||
c kinetic energy is NOT the same on exit as was on entry
|
||||
call qmd_fix_momentum(nat,dbl_mb(k_m),dbl_mb(k_v),dbl_mb(k_r),
|
||||
$ ekin,.false.)
|
||||
end if
|
||||
c
|
||||
call ga_sync()
|
||||
c
|
||||
if (ga_nodeid().eq.0) then
|
||||
final_temp = 2.0d0*ekin/kb/(1.d0*ndeg)
|
||||
write(luout,120) "Current temp. (K):",final_temp
|
||||
end if
|
||||
c
|
||||
c perform energy gradient calculation on initial coordinates
|
||||
call qmd_gradient(rtdb,nat,dbl_mb(k_g),esys)
|
||||
c
|
||||
c initialize velocity and acceleration
|
||||
call qmd_nucl_dt(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_v),dbl_mb(k_g))
|
||||
c
|
||||
c apply thermostat to initial velocities
|
||||
call qmd_thermostat(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_v),ekin)
|
||||
c
|
||||
c start nuclear molecular dynamics loop
|
||||
do istep_nucl = 1,nstep_nucl
|
||||
do istep_nucl = sstep_nucl,nstep_nucl
|
||||
c
|
||||
c update position and velocity at half step: dt2
|
||||
call qmd_nucl_dt2(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_r),
|
||||
& dbl_mb(k_v),dbl_mb(k_g))
|
||||
if (mod(istep_nucl,fstep).eq.0) then
|
||||
c remove global translations and rotations
|
||||
c kinetic energy is the same on exit as was on entry
|
||||
call qmd_fix_momentum(nat,dbl_mb(k_m),dbl_mb(k_v),dbl_mb(k_r),
|
||||
$ ekin,.true.)
|
||||
end if
|
||||
if (.not.do_none) then
|
||||
c apply thermostat to velocities
|
||||
call qmd_thermostat(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_v),ekin)
|
||||
end if
|
||||
c
|
||||
c apply thermostat
|
||||
call qmd_thermostat(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_v),ekin)
|
||||
c update velocity at half step (dt/2)
|
||||
call qmd_nucl_vel(nat,0.5d0*dt_nucl,dbl_mb(k_m),dbl_mb(k_v),
|
||||
& dbl_mb(k_g))
|
||||
c
|
||||
c update position at dt
|
||||
call qmd_nucl_pos(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_r),
|
||||
& dbl_mb(k_v))
|
||||
c
|
||||
c update geom and store in rtdb
|
||||
call qmd_geom_store(rtdb,geom,nat,dbl_mb(k_r))
|
||||
|
|
@ -149,17 +196,21 @@ c
|
|||
c calculate energy and gradient
|
||||
call qmd_gradient(rtdb,nat,dbl_mb(k_g),esys)
|
||||
c
|
||||
c update velocity at full step: dt
|
||||
call qmd_nucl_dt(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_v),
|
||||
& dbl_mb(k_g))
|
||||
c update velocity at half step (dt/2)
|
||||
call qmd_nucl_vel(nat,0.5d0*dt_nucl,dbl_mb(k_m),dbl_mb(k_v),
|
||||
& dbl_mb(k_g))
|
||||
c
|
||||
c apply thermostat
|
||||
call qmd_thermostat(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_v),ekin)
|
||||
if (.not.do_none) then
|
||||
c apply thermostat to velocities
|
||||
call qmd_thermostat(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_v),ekin)
|
||||
end if
|
||||
c
|
||||
c persist updated velocity in geom
|
||||
if (.not.geom_vel_set(geom,dbl_mb(k_v)))
|
||||
$ call errquit(pname//'geom_vel_set',0,GEOM_ERR)
|
||||
c
|
||||
c calculate current kinetic energy
|
||||
call qmd_kinetic(nat,dbl_mb(k_m),dbl_mb(k_v),ekin)
|
||||
c instantaneous total energy: instantaneous kinetic energy + system energy
|
||||
etotal = ekin + esys
|
||||
c
|
||||
|
|
@ -172,27 +223,28 @@ c
|
|||
write(luout,*)
|
||||
write(luout,*)
|
||||
call util_print_centered(LuOut,'QMD Run Information',20,.true.)
|
||||
write(luout,*) "istep_nucl: ",istep_nucl
|
||||
write(luout,*) "Time elapsed: ",istep_nucl*dt_nucl*au2fs
|
||||
write(luout,*) "ekin: ",ekin
|
||||
write(luout,*) "esys: ",esys
|
||||
write(luout,*) "etotal: ",etotal
|
||||
write(luout,*) "targ_temp: ",targ_temp
|
||||
write(luout,*) "final_temp: ",final_temp
|
||||
write(luout,*)
|
||||
write(luout,110) "Current time step:",istep_nucl
|
||||
eltime=dble(istep_nucl)*dt_nucl*au2fs
|
||||
write(luout,130) "Time elapsed (fs):",eltime
|
||||
write(luout,130) "Kinetic energy (a.u.):",ekin
|
||||
write(luout,130) "Potential energy (a.u.):",esys
|
||||
write(luout,130) "Total energy (a.u.):",etotal
|
||||
write(luout,120) "Target temp. (K):",targ_temp
|
||||
write(luout,120) "Current temp. (K):",final_temp
|
||||
|
||||
c
|
||||
c write trajectory
|
||||
call qmd_write_trj(xyzfilename,xyzfileunit,.true.,nat,
|
||||
& istep_nucl,etotal,geom)
|
||||
c
|
||||
call qmd_write_rst(rstfilename,rstfileunit,
|
||||
& nat,istep_nucl,ekin,geom)
|
||||
end if ! ga_nodeid()
|
||||
c
|
||||
end do ! istep_nucl
|
||||
c
|
||||
c write restart information
|
||||
if (ga_nodeid().eq.0) then
|
||||
call qmd_write_rst(rstfilename,rstfileunit,
|
||||
& nat,istep_nucl-1,esys,ekin,geom)
|
||||
end if
|
||||
110 format(A25,I20)
|
||||
120 format(A25,F20.2)
|
||||
130 format(A25,F20.6)
|
||||
c
|
||||
c clear memory
|
||||
if(.not.ma_pop_stack(l_act))! active
|
||||
|
|
|
|||
|
|
@ -22,6 +22,8 @@ c
|
|||
c
|
||||
logical task_gradient
|
||||
external task_gradient
|
||||
logical task_energy
|
||||
external task_energy
|
||||
c
|
||||
c preliminaries
|
||||
pname = "qmd_gradient: "
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
c
|
||||
c qmd_init: initialize parameters
|
||||
c
|
||||
subroutine qmd_init(rtdb,nstep_nucl,dt_nucl)
|
||||
subroutine qmd_init(rtdb,nstep_nucl,dt_nucl,do_none,fstep)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -22,6 +22,11 @@ c
|
|||
integer ao_bas_han ! AO basis set handle
|
||||
integer nstep_nucl ! nuclear steps
|
||||
integer nstep_elec ! total electronic steps
|
||||
integer states ! number of electronic states
|
||||
integer init_st ! initially populated state
|
||||
integer fstep ! how often to remove trans/rot
|
||||
integer tvals(8)
|
||||
integer svals
|
||||
double precision dt_nucl ! nuclear md time step
|
||||
double precision dt_elec ! electronic time step
|
||||
character*32 thermostat ! thermostat name
|
||||
|
|
@ -29,6 +34,7 @@ c
|
|||
character*32 integrator_elec ! electronic dynamics integrator
|
||||
integer ichain
|
||||
double precision nh_mass(2) ! nose-hoover mass
|
||||
logical do_none
|
||||
c
|
||||
character*32 pname
|
||||
c
|
||||
|
|
@ -45,57 +51,64 @@ c
|
|||
write(LuOut,*)
|
||||
write(LuOut,*)
|
||||
call util_print_centered
|
||||
& (LuOut,'QMD General Information',40,.true.)
|
||||
& (LuOut,'QMD Run Parameters',40,.true.)
|
||||
write(LuOut,*)
|
||||
write(LuOut,*)
|
||||
endif ! ga_nodeid
|
||||
c
|
||||
c boltzmann constant in atomic units
|
||||
kb=3.16679d-6
|
||||
kb=3.1668114d-6
|
||||
c
|
||||
c convert time in au to femtoseconds
|
||||
au2fs=0.02419d0
|
||||
au2fs=0.02418884326505d0
|
||||
c
|
||||
c convert mass in amu to au
|
||||
amu2au=1822.8884845d0
|
||||
c
|
||||
c convert length in au to Angstrom
|
||||
au2ang=0.52917721092d0
|
||||
c
|
||||
c total nuclear steps
|
||||
if (.not.rtdb_get(rtdb,'qmd:nstep_nucl',mt_int,1,nstep_nucl))
|
||||
& call errquit(pname//'failed to read nstep_nucl',0,RTDB_ERR)
|
||||
if (ga_nodeid().eq.0) write(luout,*) "nstep_nucl: ",nstep_nucl
|
||||
if (ga_nodeid().eq.0)
|
||||
& write(luout,110) "No. of nuclear steps:",nstep_nucl
|
||||
c
|
||||
c nuclear time step
|
||||
if (.not.rtdb_get(rtdb,'qmd:dt_nucl',mt_dbl,1,dt_nucl))
|
||||
& call errquit(pname//'failed to read dt_nucl',0,RTDB_ERR)
|
||||
if (ga_nodeid().eq.0) write(luout,*) "dt_nucl: ",dt_nucl
|
||||
if (ga_nodeid().eq.0)
|
||||
& write(luout,120) "Nuclear time step:",dt_nucl
|
||||
c
|
||||
c total electronic steps
|
||||
if(.not.rtdb_get(rtdb,'qmd:nstep_elec',mt_int,1,nstep_elec))
|
||||
& call errquit(pname//'failed to read nstep_elec',0,RTDB_ERR)
|
||||
if (ga_nodeid().eq.0) write(luout,*) "nstep_elec: ",nstep_elec
|
||||
c if (ga_nodeid().eq.0)
|
||||
c & write(luout,110) "No. of electronic steps:",nstep_elec
|
||||
c
|
||||
c electronic time step
|
||||
if (.not.rtdb_get(rtdb,'qmd:dt_elec',mt_dbl,1,dt_elec))
|
||||
& call errquit(pname//'failed to read dt_elec',0,RTDB_ERR)
|
||||
if (ga_nodeid().eq.0) write(luout,*) "dt_elec: ",dt_elec
|
||||
c if (ga_nodeid().eq.0)
|
||||
c & write(luout,*) "Electronic time step:",dt_elec
|
||||
c
|
||||
c target temperature
|
||||
if (.not.rtdb_get(rtdb,'qmd:targ_temp',mt_dbl,1,targ_temp))
|
||||
& call errquit(pname//'failed to read targ_temp',0,RTDB_ERR)
|
||||
if (ga_nodeid().eq.0)
|
||||
& write(luout,*) "target temperature: ",targ_temp
|
||||
c
|
||||
c nose-hoover mass parameter
|
||||
if (.not.rtdb_get(rtdb,'qmd:nh_mass',mt_dbl,2,nh_mass))
|
||||
& call errquit(pname//'failed to read nh_mass',0,RTDB_ERR)
|
||||
if (ga_nodeid().eq.0)
|
||||
& write(luout,*) "nh_mass(1): ",nh_mass(1)
|
||||
if (ga_nodeid().eq.0)
|
||||
& write(luout,*) "nh_mass(2): ",nh_mass(2)
|
||||
& write(luout,120) "Target temp. (K):",targ_temp
|
||||
c
|
||||
c thermostat
|
||||
do_berendsen = .false.
|
||||
do_nosehoover = .false.
|
||||
do_langevin = .false.
|
||||
do_BDP = .false.
|
||||
do_rescale = .false.
|
||||
do_none = .false.
|
||||
if (.not. rtdb_cget(rtdb,'qmd:thermostat',1,thermostat))
|
||||
$ call errquit(pname//'failed to read thermostat',0,RTDB_ERR)
|
||||
if (ga_nodeid().eq.0) write(luout,*) "thermostat: ",thermostat
|
||||
if (ga_nodeid().eq.0)
|
||||
$ write(luout,140) "Thermostat:",trim(thermostat)
|
||||
if (thermostat.eq.'berendsen') then
|
||||
do_berendsen = .true.
|
||||
else if (thermostat.eq.'nose-hoover') then
|
||||
|
|
@ -106,41 +119,86 @@ c thermostat
|
|||
m_nh(ichain) = nh_mass(ichain) ! mass of chains
|
||||
g_nh(ichain) = 0.d0 ! bath coupling
|
||||
end do
|
||||
else if (thermostat.eq.'langevin') then
|
||||
do_langevin = .true.
|
||||
else if (thermostat.eq.'rescale') then
|
||||
do_rescale = .true.
|
||||
else if (thermostat.eq.'BDP') then
|
||||
do_BDP = .true.
|
||||
else if (thermostat.eq.'none') then
|
||||
do_none = .true.
|
||||
else
|
||||
if (ga_nodeid().eq.0)
|
||||
& write(luout,*) "unknown thermostat"
|
||||
if (ga_nodeid().eq.0)
|
||||
& write(luout,*) "using default: berendsen thermostat"
|
||||
do_berendsen = .true.
|
||||
& write(luout,*) "using default: NO thermostat"
|
||||
do_none = .true.
|
||||
end if
|
||||
c
|
||||
c nose-hoover mass parameter
|
||||
if (.not.rtdb_get(rtdb,'qmd:nh_mass',mt_dbl,2,nh_mass))
|
||||
& call errquit(pname//'failed to read nh_mass',0,RTDB_ERR)
|
||||
if (do_nosehoover) then
|
||||
if (ga_nodeid().eq.0)
|
||||
& write(luout,130) "NH mass(1):",nh_mass(1)
|
||||
if (ga_nodeid().eq.0)
|
||||
& write(luout,130) "NH mass(2):",nh_mass(2)
|
||||
end if
|
||||
c
|
||||
c langevin friction parameter
|
||||
friction = 0.1d0
|
||||
if (.not.rtdb_get(rtdb,'qmd:friction',mt_dbl,1,friction))
|
||||
& friction = 0.1d0
|
||||
if (do_langevin) then
|
||||
if (ga_nodeid().eq.0) write(luout,130) "Friction:",friction
|
||||
end if
|
||||
c
|
||||
c berendsen tau parameter, also used for BDP
|
||||
if (.not.rtdb_get(rtdb,'qmd:tau',mt_dbl,1,tau))
|
||||
& call errquit(pname//'failed to read tau',0,RTDB_ERR)
|
||||
if (do_berendsen.or.do_BDP) then
|
||||
if (ga_nodeid().eq.0) write(luout,120) "Tau:",tau
|
||||
end if
|
||||
c
|
||||
c check dt_nucl against the berendsen tau
|
||||
if (do_berendsen.and.(dt_nucl.ge.tau)) then
|
||||
tau = 1.25d0*dt_nucl
|
||||
if (ga_nodeid().eq.0) then
|
||||
write(luout,*) "dt_nucl has to be smaller than tau"
|
||||
write(luout,*) "adjusting tau to accomodate"
|
||||
write(luout,*) "new tau: ", tau
|
||||
end if
|
||||
end if
|
||||
c
|
||||
c random number seed
|
||||
call date_and_time(values=tvals)
|
||||
c current seconds x milliseconds + minutes
|
||||
svals = tvals(7)*tvals(8)+tvals(6)
|
||||
idum = -svals
|
||||
if (.not.rtdb_get(rtdb,'qmd:rand_seed',mt_int,1,idum))
|
||||
& idum = -svals ! default
|
||||
if (ga_nodeid().eq.0) write(luout,110) "Random seed:",idum
|
||||
c
|
||||
c nuclear md integrator
|
||||
if (.not. rtdb_cget(rtdb,'qmd:integrator_nucl',1,integrator_nucl))
|
||||
$ call errquit(pname//'failed to read integrator',0,RTDB_ERR)
|
||||
if (ga_nodeid().eq.0)
|
||||
& write(luout,*) "integrator_nucl: ",integrator_nucl
|
||||
& write(luout,140) "Nuclear integrator:",trim(integrator_nucl)
|
||||
if (integrator_nucl.eq.'velocity-verlet') do_veloverlet = .true.
|
||||
c
|
||||
c berendsen tau parameter
|
||||
if (.not.rtdb_get(rtdb,'qmd:tau',mt_dbl,1,tau))
|
||||
& call errquit(pname//'failed to read tau',0,RTDB_ERR)
|
||||
if (ga_nodeid().eq.0) write(luout,*) "tau: ",tau
|
||||
c
|
||||
c check dt_nucl against the berendsen tau
|
||||
if (do_berendsen.and.(dt_nucl.ge.tau)) then
|
||||
dt_nucl = 0.75d0*tau
|
||||
if (ga_nodeid().eq.0) then
|
||||
write(luout,*) "dt_nucl has to be lesser than tau"
|
||||
write(luout,*) "adjusting dt_nucl to accomodate"
|
||||
write(luout,*) "new dt_nucl: ", dt_nucl
|
||||
end if
|
||||
end if
|
||||
c
|
||||
c is linear molecule ?
|
||||
if (.not.rtdb_get(rtdb,'qmd:linear',mt_log,1,do_linear))
|
||||
& do_linear = .false.
|
||||
c
|
||||
c remove translations and rotations every fsteps
|
||||
if (.not.rtdb_get(rtdb,'qmd:fstep',mt_int,1,fstep))
|
||||
$ fstep = 1000
|
||||
if (ga_nodeid().eq.0) call util_flush(LuOut)
|
||||
c
|
||||
110 format(A25,I20)
|
||||
120 format(A25,F20.2)
|
||||
130 format(A25,F20.6)
|
||||
140 format(A25,A20)
|
||||
c
|
||||
return
|
||||
end
|
||||
|
|
|
|||
|
|
@ -24,7 +24,10 @@ c
|
|||
integer nstep_nucl
|
||||
double precision dt_elec
|
||||
integer nstep_elec
|
||||
integer fstep
|
||||
double precision targ_temp
|
||||
double precision friction
|
||||
integer idum
|
||||
double precision tau
|
||||
double precision nh_mass(2)
|
||||
character*32 test
|
||||
|
|
@ -42,7 +45,7 @@ c Set default values
|
|||
c ------------------
|
||||
c
|
||||
c Nuclear MD time step: dt_nucl
|
||||
dt_nucl = 1.d0
|
||||
dt_nucl = 10.d0
|
||||
if (.not.rtdb_put(rtdb,'qmd:dt_nucl',mt_dbl,1,dt_nucl))
|
||||
& call errquit(pname//'failed to write dt_nucl',0,RTDB_ERR)
|
||||
c
|
||||
|
|
@ -67,12 +70,12 @@ c Target nuclear MD temperature: targ_temp (Kelvin)
|
|||
& call errquit(pname//'failed to write targ_temp',0,RTDB_ERR)
|
||||
c
|
||||
c Nuclear thermostat
|
||||
thermostat='berendsen'
|
||||
thermostat='none'
|
||||
if (.not.rtdb_cput(rtdb,'qmd:thermostat',1,thermostat))
|
||||
& call errquit(pname//'failed to write thermostat',0,RTDB_ERR)
|
||||
c
|
||||
c Berendsen tau parameter
|
||||
tau=0.1d0
|
||||
c Berendsen tau parameter, also used for BDP
|
||||
tau=1.d3
|
||||
if (.not.rtdb_put(rtdb,'qmd:tau',mt_dbl,1,tau))
|
||||
& call errquit(pname//'failed to write tau',0,RTDB_ERR)
|
||||
c
|
||||
|
|
@ -93,6 +96,11 @@ c Is linear molecule ?
|
|||
if (.not.rtdb_put(rtdb,'qmd:linear',mt_log,1,do_linear))
|
||||
& call errquit(pname//'failed to write do_linear',0,RTDB_ERR)
|
||||
c
|
||||
c Remove translations/rotations every fsteps
|
||||
fstep=1000
|
||||
if (.not.rtdb_put(rtdb,'qmd:fstep',mt_int,1,fstep))
|
||||
& call errquit(pname//'failed to write fstep',0,RTDB_ERR)
|
||||
c
|
||||
c ----------
|
||||
c Read input
|
||||
c ----------
|
||||
|
|
@ -106,8 +114,8 @@ c Nuclear MD time step: dt_nucl in femtoseconds
|
|||
if (inp_compare(.false.,test,'dt_nucl')) then
|
||||
if (.not.inp_f(dt_nucl)) then
|
||||
write(luout,*) 'dt_nucl: value not found; ',
|
||||
1 'default value of 1.0 au will be used'
|
||||
dt_nucl = 1.d0
|
||||
1 'default value of 10.0 au will be used'
|
||||
dt_nucl = 10.d0
|
||||
endif
|
||||
if (.not.rtdb_put(rtdb,'qmd:dt_nucl',mt_dbl,1,dt_nucl))
|
||||
& call errquit(pname//'failed to write dt_nucl',0,RTDB_ERR)
|
||||
|
|
@ -169,10 +177,46 @@ c Nuclear thermostat: Berendsen
|
|||
thermostat='berendsen'
|
||||
if (.not.rtdb_cput(rtdb,'qmd:thermostat',1,thermostat))
|
||||
1 call errquit(pname//'failed to write thermostat',0,RTDB_ERR)
|
||||
if (.not.inp_f(tau)) tau = 0.1d0
|
||||
if (.not.inp_f(tau)) tau = 1.d3
|
||||
if (.not.rtdb_put(rtdb,'qmd:tau',mt_dbl,1,tau))
|
||||
& call errquit(pname//'failed to write tau',0,RTDB_ERR)
|
||||
c
|
||||
c Nuclear thermostat: BDP (Bussi, Donadio, Parrinello 2007)
|
||||
else if (inp_compare(.false.,test,'BDP')) then
|
||||
thermostat='BDP'
|
||||
if (.not.rtdb_cput(rtdb,'qmd:thermostat',1,thermostat))
|
||||
1 call errquit(pname//'failed to write thermostat',0,RTDB_ERR)
|
||||
if (.not.inp_f(tau)) tau = 1.d2
|
||||
if (.not.rtdb_put(rtdb,'qmd:tau',mt_dbl,1,tau))
|
||||
& call errquit(pname//'failed to write tau',0,RTDB_ERR)
|
||||
c
|
||||
c Nuclear thermostat: Langevin (Bussi, Parrinello 2007))
|
||||
else if (inp_compare(.false.,test,'langevin')) then
|
||||
thermostat='langevin'
|
||||
if (.not.rtdb_cput(rtdb,'qmd:thermostat',1,thermostat))
|
||||
1 call errquit(pname//'failed to write thermostat',0,RTDB_ERR)
|
||||
c
|
||||
c Nuclear thermostat: Velocity rescaling
|
||||
else if (inp_compare(.false.,test,'rescale')) then
|
||||
thermostat='rescale'
|
||||
if (.not.rtdb_cput(rtdb,'qmd:thermostat',1,thermostat))
|
||||
1 call errquit(pname//'failed to write thermostat',0,RTDB_ERR)
|
||||
c
|
||||
c Friction
|
||||
else if (inp_compare(.false.,test,'friction')) then
|
||||
if (.not.inp_f(friction)) friction = 0.1d0
|
||||
if (.not.rtdb_put(rtdb,'qmd:friction',mt_dbl,1,friction))
|
||||
& call errquit(pname//'failed to write friction',0,
|
||||
& RTDB_ERR)
|
||||
c
|
||||
c Random number seed
|
||||
else if (inp_compare(.false.,test,'rand_seed')) then
|
||||
if (.not.inp_i(idum)) idum = 3989
|
||||
idum = sign(idum,-1) ! seed needs to be negative
|
||||
if (.not.rtdb_put(rtdb,'qmd:rand_seed',mt_int,1,idum))
|
||||
& call errquit(pname//'failed to write rand_seed',0,
|
||||
& RTDB_ERR)
|
||||
c
|
||||
c Nuclear MD integrator
|
||||
else if (inp_compare(.false.,test,'velocity-verlet')) then
|
||||
integrator_nucl='velocity-verlet'
|
||||
|
|
@ -185,6 +229,12 @@ c Is linear molecule ?
|
|||
if (.not.rtdb_put(rtdb,'qmd:linear',mt_log,1,do_linear))
|
||||
& call errquit(pname//'failed to write do_linear',0,RTDB_ERR)
|
||||
c
|
||||
c Remove translations and rotations every fsteps
|
||||
else if (inp_compare(.false.,test,'fstep')) then
|
||||
if (.not.inp_i(fstep)) fstep=1000
|
||||
if (.not.rtdb_put(rtdb,'qmd:fstep',mt_int,1,fstep))
|
||||
& call errquit(pname//'failed to write fstep',0,RTDB_ERR)
|
||||
c
|
||||
c Print
|
||||
else if (inp_compare(.false.,test,'print')) then
|
||||
call util_print_input(rtdb,'qmd')
|
||||
|
|
|
|||
|
|
@ -1,23 +1,25 @@
|
|||
c
|
||||
c calculate kinetic energy
|
||||
c
|
||||
subroutine qmd_kinetic(nat,m,v,ekin)
|
||||
subroutine qmd_kinetic(n,m,v,ekin)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
integer nat
|
||||
double precision m(nat)
|
||||
double precision v(3,nat)
|
||||
#include "qmd_common.fh"
|
||||
c
|
||||
integer n
|
||||
double precision m(n)
|
||||
double precision v(3,n)
|
||||
double precision ekin
|
||||
c
|
||||
integer i
|
||||
c
|
||||
c total kinetic energy
|
||||
ekin = 0.d0
|
||||
do i = 1,nat
|
||||
do i = 1,n
|
||||
ekin = ekin + m(i)*(v(1,i)*v(1,i)+v(2,i)*v(2,i)+v(3,i)*v(3,i))
|
||||
end do
|
||||
ekin = 0.5d0*ekin
|
||||
ekin = 0.5d0*ekin*amu2au
|
||||
c
|
||||
return
|
||||
end
|
||||
|
|
|
|||
|
|
@ -1,36 +1,38 @@
|
|||
c
|
||||
c qmd_nucl_dt: velocity at full dt
|
||||
c Update velocity and positions
|
||||
c
|
||||
subroutine qmd_nucl_dt(nat,dt,m,v,g)
|
||||
subroutine qmd_nucl_dt(nat,dt,m,r,v,g)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
integer nat ! number of atoms
|
||||
double precision dt ! step size
|
||||
double precision m(nat) ! mass
|
||||
double precision r(3,nat) ! coordinates
|
||||
double precision v(3,nat) ! velocities
|
||||
double precision g(3,nat) ! gradients
|
||||
c
|
||||
integer i
|
||||
double precision a(3,nat) ! acceleration
|
||||
integer i,j
|
||||
double precision dt2
|
||||
c
|
||||
character*30 pname
|
||||
c
|
||||
c preliminaries
|
||||
pname = "qmd_nucl_dt: "
|
||||
dt2 = 0.5d0*dt
|
||||
pname='qmd_nucl_dt: '
|
||||
c
|
||||
c acceleration (a) and velocities (v) at dt
|
||||
do i=1,nat
|
||||
c calculate acceleration, update positions (r) and velocities (v) at dt/2
|
||||
do i = 1,nat
|
||||
a(1,i) = -g(1,i)/m(i)
|
||||
v(1,i) = v(1,i) + a(1,i)*dt2
|
||||
v(1,i) = v(1,i) + a(1,i)*0.5d0*dt
|
||||
r(1,i) = r(1,i) + v(1,i)*dt
|
||||
c
|
||||
a(2,i) = -g(2,i)/m(i)
|
||||
v(2,i) = v(2,i) + a(2,i)*dt2
|
||||
v(2,i) = v(2,i) + a(2,i)*0.5d0*dt
|
||||
r(2,i) = r(2,i) + v(2,i)*dt
|
||||
c
|
||||
a(3,i) = -g(3,i)/m(i)
|
||||
v(3,i) = v(3,i) + a(3,i)*dt2
|
||||
v(3,i) = v(3,i) + a(3,i)*0.5d0*dt
|
||||
r(3,i) = r(3,i) + v(3,i)*dt
|
||||
end do
|
||||
c
|
||||
return
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
c
|
||||
c qmd_thermostat: apply thermostat
|
||||
c
|
||||
subroutine qmd_thermostat(nat,dt,m,v,ekin)
|
||||
subroutine qmd_thermostat(n,dt,m,v,ekin)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -12,11 +12,12 @@ c
|
|||
c
|
||||
#include "qmd_common.fh"
|
||||
c
|
||||
integer nat ! number of atoms
|
||||
double precision dt ! step size
|
||||
double precision m(nat) ! mass
|
||||
double precision v(3,nat) ! velocities
|
||||
double precision ekin ! instantaneous kinetic energy
|
||||
integer n ! number of atoms
|
||||
double precision dt ! step size
|
||||
double precision m(n) ! mass
|
||||
double precision v(3,n) ! velocities
|
||||
double precision ekin ! instantaneous kinetic energy
|
||||
double precision ekin2
|
||||
c
|
||||
character*30 pname
|
||||
c
|
||||
|
|
@ -24,13 +25,19 @@ c preliminaries
|
|||
pname='qmd_thermostat: '
|
||||
c
|
||||
c calculate kinetic energy
|
||||
call qmd_kinetic(nat,m,v,ekin)
|
||||
call qmd_kinetic(n,m,v,ekin)
|
||||
c
|
||||
c apply thermostat
|
||||
if (do_nosehoover) then
|
||||
call qmd_nosehoover(nat,dt,v,ekin)
|
||||
call qmd_nosehoover(n,dt,m,v,ekin)
|
||||
else if (do_berendsen) then
|
||||
call qmd_berendsen(nat,dt,v,ekin)
|
||||
call qmd_berendsen(n,dt,m,v,ekin)
|
||||
else if (do_langevin) then
|
||||
call qmd_langevin(n,dt,m,v,ekin)
|
||||
else if (do_rescale) then
|
||||
call qmd_rescale(n,m,v,ekin,ekin2)
|
||||
else if (do_BDP) then
|
||||
call qmd_BDP(n,dt,m,v,ekin)
|
||||
else
|
||||
call errquit(pname//'unknown thermostat',0,RTDB_ERR)
|
||||
end if
|
||||
|
|
@ -40,7 +47,7 @@ c
|
|||
c
|
||||
c qmd_nosehoover: nose-hoover thermostat
|
||||
c
|
||||
subroutine qmd_nosehoover(nat,dt,v,ekin)
|
||||
subroutine qmd_nosehoover(n,dt,m,v,ekin)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -51,9 +58,10 @@ c
|
|||
c
|
||||
#include "qmd_common.fh"
|
||||
c
|
||||
integer nat
|
||||
integer n
|
||||
double precision dt
|
||||
double precision v(3,nat)
|
||||
double precision m(n) ! mass
|
||||
double precision v(3,n)
|
||||
double precision ekin
|
||||
c
|
||||
double precision dt2,dt4,dt8
|
||||
|
|
@ -102,7 +110,7 @@ c
|
|||
v_nh(2) = v_nh(2) + g_nh(2)*dt4
|
||||
c
|
||||
c scale velocity
|
||||
do i = 1,nat
|
||||
do i = 1,n
|
||||
v(1,i) = scalefac*v(1,i)
|
||||
v(2,i) = scalefac*v(2,i)
|
||||
v(3,i) = scalefac*v(3,i)
|
||||
|
|
@ -113,7 +121,7 @@ c
|
|||
c
|
||||
c qmd_berendsen: berendsen thermostat
|
||||
c
|
||||
subroutine qmd_berendsen(nat,dt,v,ekin)
|
||||
subroutine qmd_berendsen(n,dt,m,v,ekin)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -124,9 +132,10 @@ c
|
|||
c
|
||||
#include "qmd_common.fh"
|
||||
c
|
||||
integer nat
|
||||
integer n
|
||||
double precision dt
|
||||
double precision v(3,nat)
|
||||
double precision m(n) ! mass
|
||||
double precision v(3,n)
|
||||
double precision ekin
|
||||
c
|
||||
integer i,j
|
||||
|
|
@ -142,9 +151,160 @@ c instantaneous temperature
|
|||
c
|
||||
c velocity scaling factor: Berendsen scaling
|
||||
ratio_temp = targ_temp/inst_temp
|
||||
scalefac = dsqrt(1.d0 + (dt/tau)*(ratio_temp-1.d0))
|
||||
scalefac = dsqrt(1.d0 + (dt/(2.d0*tau))*(ratio_temp-1.d0))
|
||||
c
|
||||
do i = 1,nat
|
||||
do i = 1,n
|
||||
v(1,i) = scalefac*v(1,i)
|
||||
v(2,i) = scalefac*v(2,i)
|
||||
v(3,i) = scalefac*v(3,i)
|
||||
end do
|
||||
c
|
||||
return
|
||||
end
|
||||
c
|
||||
c qmd_langevin: Langevin thermostat:
|
||||
c Bussi and Parrinello, Phys. Rev. E (2007)
|
||||
c
|
||||
subroutine qmd_langevin(n,dt,m,v,ekin)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
#include "errquit.fh"
|
||||
#include "mafdecls.fh"
|
||||
#include "stdio.fh"
|
||||
#include "global.fh"
|
||||
c
|
||||
#include "qmd_common.fh"
|
||||
c
|
||||
integer n
|
||||
double precision dt
|
||||
double precision m(n) ! mass
|
||||
double precision v(3,n)
|
||||
double precision ekin
|
||||
c
|
||||
integer i,j
|
||||
double precision inst_temp
|
||||
double precision ratio_temp
|
||||
double precision c1,c2
|
||||
c
|
||||
external qmd_gasdev
|
||||
double precision qmd_gasdev
|
||||
c
|
||||
character*32 pname
|
||||
pname='qmd_langevin: '
|
||||
c
|
||||
c calculate c1
|
||||
c1 = exp(-friction*dt/2.d0)
|
||||
do i = 1,n
|
||||
c2=sqrt((1.d0-c1**2)*(kb*targ_temp)/(m(i)*amu2au))
|
||||
do j = 1,3
|
||||
v(j,i) = c1*v(j,i)+c2*qmd_gasdev(idum)
|
||||
end do
|
||||
end do
|
||||
c
|
||||
return
|
||||
end
|
||||
c
|
||||
c
|
||||
c qmd_rescale: velocity rescaling
|
||||
c
|
||||
subroutine qmd_rescale(n,m,v,ekin,ekin2)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
#include "errquit.fh"
|
||||
#include "mafdecls.fh"
|
||||
#include "stdio.fh"
|
||||
#include "global.fh"
|
||||
c
|
||||
#include "qmd_common.fh"
|
||||
c
|
||||
integer n
|
||||
double precision m(n)
|
||||
double precision v(3,n)
|
||||
double precision ekin
|
||||
double precision ekin2
|
||||
c
|
||||
integer i
|
||||
double precision inst_temp
|
||||
double precision ratio_temp
|
||||
double precision scalefac
|
||||
c
|
||||
character*32 pname
|
||||
pname='qmd_rescale: '
|
||||
c
|
||||
c instantaneous temperature
|
||||
inst_temp = 2.0d0*ekin/kb/(1.d0*ndeg)
|
||||
c
|
||||
c velocity rescaling
|
||||
c ratio_temp = targ_temp/inst_temp
|
||||
ratio_temp = ekin/ekin2
|
||||
scalefac = dsqrt(ratio_temp)
|
||||
c
|
||||
do i = 1,n
|
||||
v(1,i) = scalefac*v(1,i)
|
||||
v(2,i) = scalefac*v(2,i)
|
||||
v(3,i) = scalefac*v(3,i)
|
||||
end do
|
||||
c
|
||||
return
|
||||
end
|
||||
c
|
||||
c
|
||||
c qmd_BDP: canonical velocity rescaling
|
||||
c Bussi, Donadio, and Parrinello
|
||||
c J. Chem. Phys. 126, 014101 (2007)
|
||||
c
|
||||
subroutine qmd_BDP(n,dt,m,v,ekin)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
#include "errquit.fh"
|
||||
#include "mafdecls.fh"
|
||||
#include "stdio.fh"
|
||||
#include "global.fh"
|
||||
c
|
||||
#include "qmd_common.fh"
|
||||
c
|
||||
integer n
|
||||
double precision dt
|
||||
double precision m(n)
|
||||
double precision v(3,n)
|
||||
double precision ekin
|
||||
c
|
||||
integer i
|
||||
double precision scheck
|
||||
double precision sigma
|
||||
double precision rr, ss, cc
|
||||
double precision scalefac
|
||||
c
|
||||
external qmd_gasdev
|
||||
double precision qmd_gasdev
|
||||
external qmd_sumno
|
||||
double precision qmd_sumno
|
||||
c
|
||||
character*32 pname
|
||||
pname='qmd_BDP: '
|
||||
c
|
||||
c idum, tau, ndeg, targ_temp in common block
|
||||
c
|
||||
if (tau.ne.0.d0) then
|
||||
cc = exp(-dt/(2.d0*tau))
|
||||
else
|
||||
cc = 0.d0
|
||||
end if
|
||||
c
|
||||
rr = qmd_gasdev(idum)
|
||||
sigma = dble(ndeg)*kb*targ_temp/2.d0
|
||||
ss = qmd_sumno(ndeg-1,idum)
|
||||
scheck = rr + dsqrt(cc*dble(ndeg)*ekin/((1.d0-cc)*sigma))
|
||||
c
|
||||
scalefac = cc + (sigma/(dble(ndeg)*ekin))*(1.d0-cc)*(ss+rr**2)
|
||||
$ + 2.d0*rr*dsqrt(cc*(1.d0-cc)*sigma/(dble(ndeg)*ekin))
|
||||
c
|
||||
scalefac = dsign(dsqrt(scalefac),scheck)
|
||||
c
|
||||
do i = 1,n
|
||||
v(1,i) = scalefac*v(1,i)
|
||||
v(2,i) = scalefac*v(2,i)
|
||||
v(3,i) = scalefac*v(3,i)
|
||||
|
|
|
|||
|
|
@ -2,7 +2,7 @@ c
|
|||
c write restart
|
||||
c
|
||||
subroutine qmd_write_rst(rstfilename,rstfileunit,
|
||||
& nat,iframe,esys,ekin,geom)
|
||||
& n,iframe,ekin,geom)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -17,14 +17,14 @@ c
|
|||
c
|
||||
character*(*) rstfilename
|
||||
integer rstfileunit
|
||||
integer nat
|
||||
integer n
|
||||
integer iframe
|
||||
double precision esys,ekin,etotal
|
||||
double precision ekin
|
||||
integer geom
|
||||
c
|
||||
integer i
|
||||
double precision r(3,nat)
|
||||
double precision v(3,nat)
|
||||
double precision r(3,n)
|
||||
double precision v(3,n)
|
||||
c
|
||||
character*32 pname
|
||||
c
|
||||
|
|
@ -35,29 +35,27 @@ c open file
|
|||
& position='rewind',status='unknown')
|
||||
c
|
||||
c write number of atoms, frame number and energy
|
||||
write(rstfileunit,110) nat
|
||||
etotal = esys + ekin
|
||||
write(rstfileunit,111) iframe,etotal,esys,ekin
|
||||
write(rstfileunit,112) targ_temp,final_temp
|
||||
write(rstfileunit,110) n
|
||||
write(rstfileunit,111) iframe,ekin
|
||||
if (.not.geom_cart_coords_get(geom,r))
|
||||
$ call errquit(pname//'geom_cart_coords_get',0,GEOM_ERR)
|
||||
if (.not.geom_vel_get(geom,v))
|
||||
$ call errquit(pname//'geom_vel_set',0,GEOM_ERR)
|
||||
c
|
||||
write(rstfileunit,*) "Coordinates"
|
||||
do i=1,nat
|
||||
do i=1,n
|
||||
write(rstfileunit,113) r(1,i),r(2,i),r(3,i)
|
||||
end do
|
||||
c
|
||||
write(rstfileunit,*) "Velocities"
|
||||
do i=1,nat
|
||||
do i=1,n
|
||||
write(rstfileunit,113) v(1,i),v(2,i),v(3,i)
|
||||
end do
|
||||
c
|
||||
if (rstfileunit.ne.0) close(unit=rstfileunit)
|
||||
c
|
||||
110 format(I5)
|
||||
111 format(I5,3f22.12)
|
||||
111 format(I10,3f22.12)
|
||||
112 format(2f22.12)
|
||||
113 format(3f22.12)
|
||||
c
|
||||
|
|
|
|||
|
|
@ -2,7 +2,7 @@ c
|
|||
c qmd_write_trj: write trajectory: xyz file in angstroms
|
||||
c
|
||||
subroutine qmd_write_trj(xyzfilename,xyzfileunit,
|
||||
& mode_append,nat,iframe,etotal,geom)
|
||||
& mode_append,n,iframe,etotal,geom)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -18,10 +18,15 @@ c
|
|||
character*(*) xyzfilename
|
||||
integer xyzfileunit
|
||||
logical mode_append
|
||||
integer nat
|
||||
integer i, j
|
||||
integer n
|
||||
integer iframe
|
||||
double precision etotal
|
||||
integer geom
|
||||
character(16) t(n)
|
||||
double precision r(3,n)
|
||||
double precision v(3,n)
|
||||
double precision q(n)
|
||||
c
|
||||
character*32 pname
|
||||
c
|
||||
|
|
@ -37,16 +42,25 @@ c open file
|
|||
end if
|
||||
c
|
||||
c write number of atoms, frame number and energy
|
||||
write(xyzfileunit,110) nat
|
||||
write(xyzfileunit,110) n
|
||||
write(xyzfileunit,111) iframe,etotal
|
||||
c
|
||||
if (.not.geom_print_xyzq(geom, xyzfileunit))
|
||||
& call errquit(pname//'geom write failed',0,GEOM_ERR)
|
||||
c get number of centers, tags, coordinates, and charges
|
||||
if (.not.geom_cart_get(geom,i,t,r,q))
|
||||
$ call errquit(pname//'geom_cart_coords_get',0,GEOM_ERR)
|
||||
c get velocities
|
||||
if (.not.geom_vel_get(geom,v))
|
||||
$ call errquit(pname//'geom_vel_set',0,GEOM_ERR)
|
||||
c
|
||||
do i=1,n
|
||||
write(xyzfileunit,4) t(i),r(:,i)*au2ang,v(:,i)
|
||||
end do
|
||||
c
|
||||
if (xyzfileunit.ne.0) close(unit=xyzfileunit)
|
||||
c
|
||||
110 format(I5)
|
||||
111 format(I5,f22.12)
|
||||
4 format(1x,a12,3f15.8,2x,3f15.8)
|
||||
110 format(I5)
|
||||
111 format(I10,f22.12)
|
||||
c
|
||||
return
|
||||
end
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue