From: Brian Radak (brian.radak_at_gmail.com)
Date: Wed Oct 03 2018 - 16:46:21 CDT
In retrospect, it might be possible in theory, but the current code almost
definitely precludes it. Flexible water is not commonly used with NAMD.
I've heard of a flexible 3 point model (maybe there are several?)but never
seen it used before.
On Wed, Oct 3, 2018, 2:14 PM Alfonso Gijón <alfonso.gijon_at_csic.es> wrote:
> 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, 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> 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