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