Re: How to simulate flexible water with TIP4P potential?

From: David Hardy (dhardy_at_ks.uiuc.edu)
Date: Wed Oct 03 2018 - 10:35:58 CDT

Brian brings up a good point about the modeling issues.

Aside from that, it happens that the code path in NAMD that repositions the lone pairs requires rigid water.

--
David J. Hardy, Ph.D.
Beckman Institute
University of Illinois at Urbana-Champaign
405 N. Mathews Ave., Urbana, IL 61801
dhardy_at_ks.uiuc.edu, http://www.ks.uiuc.edu/~dhardy/
> On Oct 3, 2018, at 9:34 AM, Brian Radak <brian.radak_at_gmail.com> wrote:
> 
> Hi Alfonso,
> 
> I don't think that NAMD can do this, nor am I sure that it makes sense. How can you place a bisector lone pair on an asymmetric water molecule? The rigid constraints and lone pairs should be separate parts of the code, but there may be specific assumptions that a flexible model violates.
> 
> HTH, 
> BKR
> 
> 
> On Wed, Oct 3, 2018, 9:16 AM Alfonso Gijón <alfonso.gijon_at_csic.es <mailto:alfonso.gijon_at_csic.es>> wrote:
> 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 Nov 18 2019 - 23:20:08 CST