From: Alfonso Gijón (alfonso.gijon_at_csic.es)
Date: Wed Oct 03 2018 - 08:15:12 CDT
Hi, I am trying to simulate flexible water with TIP4P potential. I have
already simulated rigid water with TIP4P successfully, and I have got an
equilibrated system at 350 K. From these equilibrated positions and
velocities, I try to start a new simulation just writing "rigidBonds
none" instead of "rigidBonds all" in the configuration file of NAMD, but
energy is not conserved anymore, and after a few timesteps the program
fails.
Even if I use a timestep smaller than 1 fs (0.1 or 0.01 fs) energy is
not conserved. Any idea about what is happening and how I could simulate
flexible water? Should I generate a new structure file (.psf) for
flexible water or should I fix the distance between the lonepair atom
(whose mass is zero) and the oxigen atom? Below this message, I attached
the configuration file that I am using for my simulations.
Thanks in advance,
Alfonso
#############################################################
## JOB DESCRIPTION
#############################################################
#
# NVE simulation of water molecules
#
#############################################################
## ADJUSTABLE PARAMETERS
#############################################################
#set inputname
set outputname nve
set restart 0 ;#0 for new sim, 1 for continuation
proc get_first_ts { xscfile } {
set fd [open $xscfile r]
gets $fd
gets $fd
gets $fd line
set ts [lindex $line 0]
close $fd
return $ts
}
#set temperature 350
structure positions_autopsf.psf
coordinates equilibrated-rigid.coor
velocities equilibrated-rigid.vel
if {$restart == 1} {
bincoordinates $inputname.restart.coor
binvelocities $inputname.restart.vel
extendedSystem $inputname.restart.xsc
set currenttimestep [get_first_ts $inputname.restart.xsc]
} else {
#temperature $temperature
set currenttimestep 0
}
firsttimestep $currenttimestep
#############################################################
## SIMULATION PARAMETERS
#############################################################
# Input
paraTypeCharmm on
parameters tip4p.par
waterModel tip4
# Force-Field Parameters
exclude scaled1-4
1-4scaling 1.0
cutoff 10.0
#cutoff 15.0
switching on
switchdist 8.0
pairlistdist 12
#pairlistdist 20
# Integrator Parameters
timestep 1.0
nonbondedFreq 1
fullElectFrequency 1
stepspercycle 10
# Rigid bonds
if {1} {
# rigidBonds all
rigidBonds none
}
#Constraints and restraints
if {0} {
constraints on
consref file
conskfile file
conskcol B
}
if {0} {
fixedAtoms on
fixedAtomsFile file
fixedAtomsCol O
}
# Lowe-Andersen thermostat
#loweAndersen on
#loweAndersenTemp 300
#loweAndersenRate 100
# rescaling velocities
#rescaleFreq 1
#rescaleTemp 350
# Periodic Boundary Conditions
if {$restart == 0} {
cellBasisVector1 124 0 0
cellBasisVector2 0 124 0
cellBasisVector3 0 0 124
cellOrigin 0 0 0
}
#PME electrostatics
PME yes
PMEGridSizeX 124
PMEGridSizeY 124
PMEGridSizeZ 124
wrapAll off
# Constant Pressure Control (variable volume)
useGroupPressure yes ;# needed for rigidBonds
useFlexibleCell no
useConstantArea no
# Output
outputName $outputname
if {1} {
binaryoutput no
#restartfreq 1000
#dcdfreq 1000
#xstFreq 1000
outputEnergies 1
outputPressure 1
} else {
binaryoutput no
#restartfreq 1000
#dcdfreq 1000
#xstFreq 1000
outputEnergies 1
outputPressure 1
}
#############################################################
## EXECUTION SCRIPT
#############################################################
# Minimization
#if {$restart == 0} {
#minimize 500
#reinitvels $temperature
#}
run 300
This archive was generated by hypermail 2.1.6 : Mon Dec 31 2018 - 23:21:25 CST