Re: tclBC during SMD and ABF simulations

From: Jérôme Hénin (jerome.henin_at_ibpc.fr)
Date: Mon Nov 28 2016 - 12:45:48 CST

Hi Olga,

Right now you seem to be pushing ions away from a fixed position. That
should be fine as long as that position is correct, that is unless the
protein drifts away.

If you were to use the position of the tagged ion instead, and apply ABF to
that ion, then things would get more complicated because your restraint
would introduce a coupling between the tagged ions and the others.

Jerome

On 28 November 2016 at 19:36, Olya Kravchenko <ovkrav_at_gmail.com> wrote:

> Hi all,
>
> I would appreciate if someone could comment on the script I made to push
> ions away from the tagged ion during abf simulations inside the ion
> channel. It seems to work the way it is now for simple equilibration, but I
> am not sure whether this is correct, especially when combined with SMD
> and/or ABF. If I apply force to the tagged ion, can it interfere somehow
> with the tclBC that I apply to surrounding ions?
>
> The simulation system is a protein channel solvated in NaCl and the
> permeating (tagged) ion is sodium with a known atom id number (344742,
> obtained by get serial command in VMD). I want to push all ions away from
> the tagged ion and leave the protein unaffected. For that I distinguish
> ions from the protein's atoms by mass. The sphere radius and the force
> strength are set in the main configuration file. Here is a tclBC scrpt that
> is called from within the main configuration file:
>
> #Set a list of atom IDs of all the ions
>
> set ionList {}
> set atomID 1
> set inStream [open $pdbSource r]
> set sphereCenter {22.5 89.6 93.2} ; #this is the position of the tagged
> ion
>
> # go through pdb file and find entries for ions and save them to the list
>
> foreach line [split [read $inStream] \n] {
>
> set string1 [string range $line 0 3]
> set string2 [string range $line 17 20]
> set string2 [string trim $string2]
>
> if { [string equal $string1 {ATOM}] || \
> [string equal $string1 {HETA}] } {
>
> if { [string equal $string2 {CLA}] || \
> [string equal $string2 {SOD}] } {
> lappend ionList $atomID
> }
>
> incr atomID
> }
> }
> close $inStream
>
> wrapmode cell
>
> proc calcforces {step unique} {
>
> global atomID sphereCenter sphereRadius maxForce ionList
>
> # find the components of the sphere's center
>
> foreach { x0 y0 z0 } $sphereCenter { break }
>
> while {[nextatom]} {
>
> set atomid [getid] ;# get the ID of the current atom
>
> # check if this ID is listed in ionList
>
> if { [lsearch $ionList $atomid] >= 0 } {
>
> set rvec [getcoord] ;# get ion's coordinates
> foreach { x y z } $rvec { break }
>
> # get distance between the ion and the sphere's center
>
> set rho [expr sqrt(($x-$x0)*($x-$x0) + ($y-$y0)*($y-$y0) + \
> ($z-$z0)*($z-$z0))]
>
> if { $atomid == 344742 || [getmass] < 22 || [getmass] > 36 } {
> dropatom
> continue
>
> } else {
>
> set forceX [expr $maxForce * ($x-$x0) / $rho]
> set forceY [expr $maxForce * ($y-$y0) / $rho]
> set forceZ [expr $maxForce * ($z-$z0) / $rho]
> }
>
> addforce "$forceX $forceY $forceZ"
> }
> }
> }
>
> Does it look correct? I would very much appreciate your comments! Thanks
> in advance,
>
> Olga
>

This archive was generated by hypermail 2.1.6 : Sun Dec 31 2017 - 23:20:50 CST