Re: How to simulate flexible water with TIP4P potential?

From: Alfonso Gijón (alfonso.gijon_at_csic.es)
Date: Wed Oct 03 2018 - 13:13:59 CDT

Thanks for the answers. I think that it would be always possible to
place the massless M atom along the symmetry axis between the hydrogen
and oxygen atoms, maintaining the distance OM constant, even when
stretching (for OH distance) and bending (for HOH angle) motion is
included. But maybe NAMD is not able to work with these features.

Anyway, do you know any other flexible model for water which works on NAMD?

Thanks,

Alfonso

On 03/10/18 17:35, David Hardy wrote:
> 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 <mailto:dhardy_at_ks.uiuc.edu>,
> http://www.ks.uiuc.edu/~dhardy/ <http://www.ks.uiuc.edu/%7Edhardy/>
>
>> On Oct 3, 2018, at 9:34 AM, Brian Radak <brian.radak_at_gmail.com
>> <mailto: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 : Tue Dec 31 2019 - 23:20:14 CST