Flexible Water with Multiple Time Steps and Beadiabatic:
--------------------------------------------------------
dvdr holds low frequency forces
dvdr2 holds high frequency forces.
RP Contraction:
----------------
Splits the LJ and ewald from the intramolecular force and uses RP contraction on them.
07-03-08: Split LJ into a new subroutine. Beadiabatic forces routine is now passed iopt which determines which part of the force routine to perform beadiabatic separation on.
3 Split (3s) version : uses a different number of normal modes for the ewald sum as for the LJ.
10_03_07 : added rg routine static properties calculation
10_03_07_05 : beadd stores and repaces r rather than backwards FT.
11_03: added normal mode as primitive represenation (evolve_nm) which dramatically speeds up the beadiabatic program when used with large multiple time steps (e.g. MTS=5 or MTS=10) since it reduces the FTs from 44 per timestep to 23 at MTS=1 and 8 per timestep to 5 per timestep at MTS = 1. Beadiabatic evolution also incorporated into normal mode primtive scheme using beadd_forces_nm subroutine which takes in positions in normal mode rep and produces normal mode forces.
13_03 : RDF speeded up
13_03_p : Program now reads parameters from a parameter file. Code tidied up using ftnchek. mts is now in input file.
16_03 : Output of <vinter>,<vintra>,<tintra>,<tinter> seperately. Improved handling of exact estimators.
02_04: Added output of ewald and LJ energies seperately.
03_04: CHANGE IN ORDERING OF BEADD NORMAL MODE FREQUENCIES
14_05 : Program code cleaned up.
20_04: OCF for Classical fixed. ACMD evolution added.
01_05: Smeared Charge Ewald intermolecular and Fourth Order Morse intramolecular added.
06-05-f : Ice starting configuration. FCC initial lattice and non-cubic ewald.
12-05 : OCF problems fixed
30-05 : Linked cell list for LJ - used when nm = 700 or above.
09-06 : Major speed improvements in LJ evaluation and electrostatic contraction scheme. For large systems a linked cell list is used to create the neighbour list.
26-06 : Added knowledge of M-site to electrostatic RP contraction. Note:// still cannot use smeared charges with electrostatic contraction scheme in this version. Renamed to water_'date' from water_bea_nm_'date'-demwald.
03-07 : Minor changes. Param file now contains masses.
08-07 : tem RP electrostatic contraction implemented. Replacing dem method (although routines still included)
09-07 : Minor fixes to problems that stopped earlier versions comiling correctly using g77. Final version written in g77 :-(.
09-07-f90 : Program converted to take advantage of fortran 90. Allocatable arrays used allowing the use of large system sizes. Globals.inc no longer contains fixed size array parameters.
3-08 : Non-cubic linked cell list for non-cubic LJ and rwald.
4-08 : Non-cubic smeared charges added. Pressure calculation and Berendsen barostat.
5-08 : Cubic pressure calculation added. RP contraction NPT now working.
31-08 : Interface setup, Non-isotropic berendsen barostat, Option for harmonic bend with anharmonic stretch.
05-09 : Rigid water evolution added for NVT. Rigid virial not yet complete.
08-09 : Rigid NPT complete.
09-09 : Buckingham O-O interaction added.
16_09im : A few small problems removed and unneccesary variables culled.
13-10 : Dielectric constant calculation error corrected. RDF calculation and exact estimators now controlled by switches in input file.
14-10 : Improved SHAKE/RATTLE routine. Also fixed crashingof classical simualtions when beadiabatic paramaters set and density profile error in the melting routine.
16-10 : Parinello global thermostat implemented for classical rigid and flexible simulations. Parinello local thermostat also availiable : for relative merits see Comput. Phys. Comm. 179 (2008) p26-29.
5-11 : Additional NPT and pressure properties calculated by md_static routine. Fixed bug in rigid classical barostat.
29-01-09 : MAJOR sampling fix done!
23-02-09 : Changed to Fortran 90
Possible Future additions:
---------------------------
1.) Inclusion of Solute atoms/molecule.
2.) RP electrostatic contraction for smeared charges.
3.) NPT simulation / pressure calculation.
4.) Polarizability.
5.) Generalize solvent treatment.
Input parameters:
------------------
ens - ensemble - option of NVT or NPT
isotope - Type of molecule (D2O, HDO, H2O (Default))
temp - temperature in K
pres - pressure in bar
rho - density in g cm**(-3)
lattice - Starting lattice : CUB, FCC, ICE or VAC
vacfac - Factor for Surface (VACuum) Calculation ncellxyz = x y z*vacfac
iamcub - Allows selection of non-cubic vs cubic energy evaluation
dtfs - time step in fs
ecut - Ewald sum accuracy parameter. Set to 0.d0 for automatic
selection of ewald parameters
nt - Dynamic simulation trajectory length
ne - Equilibration steps
ntherm - Number of Steps between parallel trajectories
nb - Number of beads
m - Number of dynamic trajectories
ng - Number of static property sampling steps
print - Integers controlling which properties should be printed/saved (4 integers, 1 or 0: Positions, Forces, Velocities, molecular dipolemoment)
pt - Print trajectory interval
pb - Print beads interval
ncellxyz - Number of cells in each direction for initialization (3 integers)
irun - Random number seed (integer)
itcf - Integers controlling calculation of dynamic properties (3 integers, 1 or 0: Cvv, Dipole Spectrum (IR), OCF)
itst - Integers controlling calculation of static properties (2 integers, 1 or 0: RDF, Exact Estimators)
rcut - Real-space cut-off for LJ interactions in Angstroms
type - Type of calculation: RPMD or ACMD
therm - Thermostat - NON(None),AND(Anderson),PRL(Parrinello
Local),PRG(Parrinello Global), (....also GLE!)
ttaufs - tau parameter for Parrinello thermostats
baro - Barostat, 'BER' or 'MCI'
taufs - tau parameter for barostat
mts - Number of multiple time steps
om - ACMD adiabatic separation parameter.
nbdf1 - RP contraction modes for ewald forces (set to zero for full nb, set to 1 for electrostatic contraction)
nbdf2 - RP contraction modes for LJ forces (set to zero for full nb)
nbdf3 - RP contraction modes for RPMDDFT forces
sig - Electrostatic RP Contraction cut-off (Angstroms, usually ~ 5.0)
rpmddft - Use RPMDDFT force Calculation
rpmde3b - ALso calculate E3B (in addition to the usual force field). Does not work with RPMDDFT.
rctdk - Use the T.D. Kühne Method (important nbdf1 and nbdf2 have to be 0, nbdf3 is the delta Force contraction)
aieq - Use RPMDDFT forces for equilibration too
Parameters for popular water model
-----------------------------------
SPC/F
------
¶m
wmass = 32831.2525d0
omass = 29156.9471d0
hmass = 1837.1527d0
qo = -0.82d0
alpha = 1.d0
oo_sig = 5.98098d0
oo_eps = 2.477083d-4
oo_gam = 0.d0
thetad = 109.47d0
reoh = 1.889726d0
de = 0.d0
alp = 0.d0
deb = 0.d0
alpb = 0.d0
wm = 0.d0
wh = 0.d0
&end
SPCFW
------
¶m
wmass = 32831.2525d0
omass = 29156.9471d0
hmass = 1837.1527d0
qo = -0.82d0
alpha = 1.d0
oo_sig = 5.981913d0
oo_eps = 2.47686d-4
oo_gam = 0.d0
thetad = 113.24d0
reoh = 1.912403d0
apot = 0.472655d0
bpot = 0.120954d0
alp = 0.d0
alpb = 0.d0
wm = 0.d0
wh = 0.d0
&end
TIP3P/FS
--------
¶m
wmass = 32831.2525d0
omass = 29156.9471d0
hmass = 1837.1527d0
qo = -0.834d0
alpha = 1.d0
oo_sig = 5.95377d0
oo_eps = 2.42542d-4
oo_gam = 0.d0
thetad = 104.5d0
reoh = 1.814137d0
apot = 0.4726468d0
bpot = 0.10850d0
alp = 0.d0
alpb = 0.d0
wm = 0.d0
wh = 0.d0
&end
QSPCFW
-------
¶m
wmass = 32831.2525d0
omass = 29156.9471d0
hmass = 1837.1527d0
qo = -0.84d0
alpha = 1.d0
oo_sig = 5.981913d0
oo_eps = 2.47686d-4
oo_gam = 0.d0
thetad = 112.d0
reoh = 1.889726d0
apot = 0.472655d0
bpot = 0.120954d0
alp = 0.d0
alpb = 0.d0
wm = 0.d0
wh = 0.d0
&end
Scott's potential
------------------
¶m
wmass = 32831.2525d0
omass = 29156.9471d0
hmass = 1837.1527d0
qo = -1.084d0
alpha = 0.797d0
oo_sig = 5.761d0
oo_eps = 3.5211d-4
oo_gam = 0.d0
thetad = 106.416d0
reoh = 1.77956d0
apot = 0.17d0
bpot = 0.52d0
alp = 1.22d0
alpb = 0.37d0
wm = 0.482d0
wh = 0.686d0
&end
TIP4P/05 rigid
---------------
¶m
wmass = 32831.2525d0
omass = 29156.9471d0
hmass = 1837.1527d0
qo = -1.1128d0
alpha = 0.73612d0
sigma = 5.96946d0
epslon = 2.95147d-4
thetad = 104.52d0
reoh = 1.8088465d0
apot = 0.d0
bpot = 0.d0
alp = 0.d0
alpb = 0.d0
wm = 0.d0
wh = 0.d0
&end
TIP4P/05 flexible - Scott
--------------------------
¶m
wmass = 32831.2525d0
omass = 29156.9471d0
hmass = 1837.1527d0
qo = -1.1128d0
alpha = 0.73612d0
oo_sig = 5.96946d0
oo_eps = 2.95147d-4
oo_gam = 0.d0
thetad = 107.4d0
reoh = 1.78d0
apot = 0.185d0
bpot = 0.07d0
alp = 1.21d0
alpb = 0.d0
wm = 0.d0
wh = 0.d0
&end
LJ :
-----
TIP4P/2005 :
epslon = 2.95147d-4
sigma = 5.96946d0
TIP4P :
epslon = 2.47012d-4
sigma = 5.96020d0
SPC :
epslon = 2.477083d-4
sigma = 5.98098d0
q-SPC/Fw :
epslon = 2.4768586d-4
sigma = 5.981913d0
Charges and alpha :
--------------------
TIP4P/2005 :
qo = -1.1128d0
alpha = 0.73612d0
TIP4P :
qo = -1.04d0
alpha = 0.74398d0
SPC/F :
qo = -0.82d0
alpha = 1.d0
q-SPC/Fw :
qo = -0.84d0
alpha = 1.d0