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 : Mon Dec 31 2018 - 23:21:26 CST