Re: question regarding simulated annealing

From: Chris Harrison (charris5_at_gmail.com)
Date: Mon Sep 10 2012 - 23:17:47 CDT

Puspita,

Puspita Halder <puspitah_at_gmail.com> writes:
> Hi Chris,
>
> Thanks for your mail. I've heard about the term vacuum bubble , is this
> totally devoid of any solvent molecule?

Yes

> I guess my situations is not
> exactly so as I already mentioned that I've got a non-uniform water
> distribution around the protein in the solvation box. Sorry, I didn't
> measure the density of water molecules in the different region of the box
> yet. I used some tcl scripts earlier but I am not so familiar with writing
> such scripts. Do you have any script available what I can use for this
> type of density calculation?

See below my sig. Put the script in a text file. Load in to the tkConsole with "source filename.tcl". Call with: checkWatCount -sel $mysel -xgridsize numberCellsInX -ygridsize numCellsInY -zgridsize numCellsInZ

If the hflag is set to 0 (ie -hflag 0), which is default, then hydrogens are not included in the operation (which seems more sensible to me). Set to "1" includes the hydrogens in the operation.

If the -output flag is set then you can dump the output to a filename, ie -output myOutputFile

All flags are optional. Default of -sel is "all"; for -xgridsize, -ygridsize, and -zgridsize it is "1", for -hflag it is "0" which indicates to not include hydrogens, and for -output the default is to the tkConsole.
 
Best,
Chris

proc checkWatCount { args } {

  foreach {i j} $args {
    if {$i=="-sel"} then { set sel $j }
    if {$i=="-xgridsize"} then { set numgridx $j }
    if {$i=="-ygridsize"} then { set numgridy $j }
    if {$i=="-zgridsize"} then { set numgridz $j }
    if {$i=="-hflag" && $j > 0} then { set hvar ""}
    if {$i=="-output"} then { set otarg $j }
  }
  if { ![info exists sel] } then {
    set sel [atomselect top all]
  }
  if { ![info exists numgridx] } then { set numgridx 1 }
  if { ![info exists numgridy] } then { set numgridy 1 }
  if { ![info exists numgridz] } then { set numgridz 1 }
  if { ![info exists hvar] } then { set hvar "noh and" }
  if { ![info exists otarg] } then {
    set ostm "stdout"
  } else {
    set ostm [open $otarg w[
  }
  
  set watMass 18.02 ; #g/m
  set watMassKg .01802 ; #kg/m
  set NA [expr {6.02214078*[expr {pow(10,23)}]}]; #Avogadro
  set a2m [expr {pow(10,-10)}]
  set a2m3 [expr {pow($a2m,3)}]

  set mymolid [$sel molid]
  set minmax [measure minmax $sel]
  set vec [vecsub [lindex $minmax 1] [lindex $minmax 0]]
  set cellxsize [lindex $vec 0]
  set cellysize [lindex $vec 1]
  set cellzsize [lindex $vec 2]
  set gridspacex [expr $cellxsize / $numgridx]
  set gridspacey [expr $cellysize / $numgridy]
  set gridspacez [expr $cellzsize / $numgridz]

  puts $ostm [format "Cell Size --> x=%4.4g y=%4.4g z=%4.4g" $cellxsize $cellysize $cellzsize]
  puts $ostm [format "Grid Spacing --> x=%4.4g y=%4.4g z=%4.4g" $gridspacex $gridspacey $gridspacez]
  puts $ostm ""
  puts $ostm "Cell\t#Waters\t%Water\tnumberDensity\tAmagat\tmolarConcentration\tdensity"
  puts $ostm "\t\t\tm^-3\t\t\tmol/(m^3)\t\tkg/(m^3)"

  set xmin [lindex [lindex $minmax 0] 0]
  for { set xind 1 } { $xind <= $numgridx } {incr xind} {
    set xmax [expr $xmin + $gridspacex ]
    set ymin [lindex [lindex $minmax 0] 1]
    for { set yind 1 } { $yind <= $numgridy } {incr yind} {
      set ymax [expr $ymin + $gridspacey ]
      set zmin [lindex [lindex $minmax 0] 2]
      for { set zind 1 } { $zind <= $numgridz } {incr zind} {
        set zmax [expr $zmin + $gridspacez ]
        set selcellall [atomselect top "$hvar x>$xmin and x<$xmax and y>$ymin and y<$ymax and z>$zmin and z<$zmax"]
        set selcellwat [atomselect top "$hvar water and element O and x>$xmin and x<$xmax and y>$ymin and y<$ymax and z>$zmin and z<$zmax"]
        
        set watcnt [$selcellwat num]
        set allcnt [$selcellall num]
        if {$watcnt == 0} {
          puts $ostm "Cell $xind $yind $zind contains 0 waters."
          set zmin [expr $zmin + $gridspacez ]
          break
        }
        if {$allcnt == 0} {
          puts $ostm "Cell $xind $yind $zind contains 0 atoms."
          set zmin [expr $zmin + $gridspacez ]
          break
        }
        set volAng [expr ($xmax - $xmin) * ($ymax - $ymin) * ($zmax - $zmin) ]
        set volMeter [expr ($xmax - $xmin) * ($ymax - $ymin) * ($zmax - $zmin) * $a2m3 ]
        set numDenAng [expr $watcnt / $volAng]
        set numDenMeter [expr $watcnt / $volMeter ]
        set numDenAmagat [expr $numDenMeter / (2.6867774 * [expr {pow(10,25)}])]
        set molarConc [expr {($numDenMeter / $NA)} ]
        set dens [expr {(($watcnt * $watMass) / ($NA * 1000)) / ($volMeter)}]
        set watpercent [expr ($watcnt / $allcnt)*100 ]
        puts $ostm [format "$xind $yind $zind\t$watcnt\t$watpercent\t%2.4g\t\t%4.4g\t%2.4g\t\t%0.4f" $numDenMeter $numDenAmagat $molarConc $dens]
        set zmin [expr $zmin + $gridspacez ]
      }
      set ymin [expr $ymin + $gridspacey ]
    }
    set xmin [expr $xmin + $gridspacex ]
  }
}

>
> Thanks for your help.
>
> Puspita
>
>
>
>
> On Thu, Aug 30, 2012 at 6:45 PM, Puspita Halder <puspitah_at_gmail.com>
>
> >
> > Hi Norman,
> >
> > Thanks for your mail. So your suggestion is to carry out the simulation at
> > 300K under NPT condition since I started my run from one of the structures
> > taken from a 500K trajectroty under the same condition . I am little
> > confused here. Did you mean to say that running the simulation at 300K
> > under NPT followed by annealing under NVT or running that under NPT as a
> > whole?
> >
> > I didn't get the point of different orientated water layers or hydrophobic
> > and hydrophilic effects. Can you clarify this a little more? I didn't
> > measure the water density exactly but by loading the trajectory in vmd it
> > seemed that water distribution around the protein is not uniform.
> >
> > Thanks again.
> >
> > Puspita
> >
> >
> >
> > On Thu, Aug 30, 2012 at 1:26 PM, Norman Geist <
> > norman.geist_at_uni-greifswald.de> wrote:
> >
> >> Also, ****
> >>
> >> ** **
> >>
> >> if the method you measure the density with, is really correct, couldn’t
> >> it be that your proteins hydrate shell is just like that because of
> >> different orientated water layers or hydrophobic and hydrophilic effects?
> >> ****
> >>
> >> ** **
> >>
> >> Norman Geist.****
> >>
> >> ** **
> >>
> >> *Von:* Norman Geist [mailto:norman.geist_at_uni-greifswald.de]
> >> *Gesendet:* Donnerstag, 30. August 2012 08:50
> >> *An:* 'Puspita Halder'
> >> *Cc:* Namd Mailing List (namd-l_at_ks.uiuc.edu)
> >> *Betreff:* AW: namd-l: question regarding simulated annealing****
> >>
> >> ** **
> >>
> >> Hi,****
> >>
> >> ** **
> >>
> >> I could imagine that the NPT run with a temperature of 500K generates
> >> slightly other box dimensions as with 300K. So maybe you would have to NPT
> >> again under the 300K conditions.****
> >>
> >> ** **
> >>
> >> Norman Geist.****
> >>
> >> ** **
> >>
> >> *Von:* owner-namd-l_at_ks.uiuc.edu [mailto:owner-namd-l_at_ks.uiuc.edu] *Im
> >> Auftrag von *Puspita Halder
> >> *Gesendet:* Dienstag, 28. August 2012 14:04
> >> *An:* namd-l_at_ks.uiuc.edu
> >> *Betreff:* namd-l: question regarding simulated annealing****
> >>
> >> ** **
> >>
> >> Dear NAMD users,
> >>
> >>
> >> I have a question regarding the water density in the solvation box during
> >> a simulated annealing run. I have started the simulated annealing run of my
> >> protein from a structure taken from a high temperature (500K) trajectory of
> >> the same generated after the minimization, heating and md run for ~2 ns at
> >> 500K temperature under NPT condiition. I decreased the temperature of my
> >> system staring from 500 K to 300 K in 5K temperature steps and for each and
> >> individual temperature step I equilibrated the system for 100 ps and
> >> finally equilibrated the system at 300K for 2 ns under NVT condition. The
> >> annealing run was seemed to be ok from the log file . The only problem that
> >> I found is the water density around the protein in the solvation box. This
> >> is found to be kind of inhomogeneous means in some portion of the box water
> >> is more dense whereas in the other portion it is less dense though the
> >> protein remains always solvated . Is this due to the effect of temperature
> >> decrease or something else? I have checked the xst file for the dimension
> >> of the box and it was fixed.
> >>
> >> If you need my config file for the simulated annealing run I'd send that

Best,
Chris

--
Chris Harrison, Ph.D.
NIH Center for Macromolecular Modeling and Bioinformatics
Theoretical and Computational Biophysics Group
Beckman Institute for Advanced Science and Technology
University of Illinois, 405 N. Mathews Ave., Urbana, IL 61801
http://www.ks.uiuc.edu/Research/namd       Voice: 773-570-6078
http://www.ks.uiuc.edu/~char               Fax:   217-244-6078

This archive was generated by hypermail 2.1.6 : Tue Dec 31 2013 - 23:22:33 CST