How to simulate flexible water with TIP4P potential?

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 Sep 16 2019 - 23:20:01 CDT