From: Azadeh Alavizargar (azadeh.alavizargar_at_gmail.com)
Date: Fri Jan 13 2017 - 15:00:16 CST
Hello Olya
Thank you very much for your explanations and the code. It is definitely
helpful. May I know what "maxforce" did you use? Actually, before I send
the email, I had prepared a code for this purpose but in the form of
TCLforces rather than TclBC, which I have copied below. I have considered a
cylinder, instead of sphere, whose axis passes through the center of mass
of some alpha carbons of the pore. It is straightforward to change it to a
sphere. When I run this code, it does not seem to slow down the simulations
much. I am not completely sure though.
I hope it is helpful anyway.
Best,
Azadeh
# atoms selected for calculating center of mass of the some alpha carbons
of the pore
set ri1 554
set ri2 3646
set ri3 6738
set ri4 9830
set grp1 {}
lappend grp1 $ri1
lappend grp1 $ri2
lappend grp1 $ri3
lappend grp1 $ri4
foreach atom $grp1 {
addatom $atom
}
# selecting ions except the target ion
set grp2 {}
for { set ii 131668 } { $ii <= 131730 } { incr ii } {
lappend grp2 $ii
}
for { set jj 131732 } { $jj <= 131754 } { incr jj } {
lappend grp2 $jj
}
foreach atom $grp2 {
addatom $atom
}
set s [addgroup $grp2]
# parameters for the cylinder and the restraining force
set cyl_r 10.0 # the cylinder radius
set cyl_k 100.0 # the force constant
proc calcforces {} {
global grp1
global grp2
global cyl_r
global cyl_k
# get coordinates and masses
loadcoords coordinate
loadmasses masses
# calculating center of mass of some alpha carbons
set s_sum "0.0 0.0 0.0"
set totalmass 0.0
foreach atom $grp1 {
set tmp [vecscale $masses($atom) $coordinate($atom)]
set s_sum [vecadd $s_sum $tmp]
set totalmass [expr $totalmass + $masses($atom)]
}
set invmass [expr 1.0/$totalmass]
set com [vecscale $invmass $s_sum]
set com_x [lindex $com 0]
set com_y [lindex $com 1]
set com_z [lindex $com 2]
# restraining ions through a cylinder
foreach atoms $grp2 {
set st $coordinate($atoms)
set sx [lindex $st 0]
set sy [lindex $st 1]
set sz [lindex $st 2]
if { $sz > 3 } {
set r2 [expr ($sx - $com_x) * ($sx - $com_x) + ($sy - $com_y) *
($sy - $com_y)]
if { $r2 < ($cyl_r*$cyl_r) } {
set sr [expr sqrt($r2)]
set fs [expr -2.*$cyl_k*($sr-$cyl_r)/$sr]
set fsx [expr $sx*$fs/$sr ]
set fsy [expr $sy*$fs/$sr ]
set fsz 0
addforce $atoms "$fsx $fsy $fsz"
}
}
}
}
On Thu, Jan 12, 2017 at 8:09 PM, Olya Kravchenko <ovkrav_at_gmail.com> wrote:
> Hi Azadeh,
>
> I recently asked similar question here about restricting the flow of ions
> through the pore in a convex surface of a spherical protein, you can try to
> find this discussion. I ended up setting up an "exclusion sphere" with a
> centre near the middle of the channel. The ions that enter the sphere are
> pushed out except for the one I want to stay inside.
>
> I used the "User-Defined Forces in NAMD" tutorial: http://www.ks.uiuc.edu/
> Training/Tutorials/science/forces/forces-tutorial-html/
> forces-tutorial.html
> particularly this section: http://www.ks.uiuc.edu/
> Training/Tutorials/science/forces/forces-tutorial-html/node3.html
> and I based my script on sample scripts that come along with it:
> http://www.ks.uiuc.edu/Training/Tutorials/science/
> forces/forces-tutorial-files/tclBCfiles/
>
> I noticed that it significantly slows down my calculations. If you come up
> with a more efficient way of removing the ions, please share.
>
> Here is the script that works in my system:
>
> set ionList {} ;# a list of atom IDs of the ions
> set atomID 1
>
> set inStream [open $pdbSource r] ;# open PDB file
>
> #set sphereCenter {34.0 98.4 92.3} ;# this sets can be set in the
> configuration file
>
> # find all sodium and chloride ions in the pdb file and add their IDs into
> 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
>
> foreach { x0 y0 z0 } $sphereCenter { break }
>
> while {[nextatom]} {
>
> set atomid [getid]
>
> # check whether atom ID is on the list, then check it's position
>
> if { [lsearch $ionList $atomid] >= 0 } {
>
> set rvec [getcoord]
> foreach { x y z } $rvec { break }
>
> # find the 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 { [getid] == 344742 || [getid] < 344526 } {
> 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"
> }
> }
> }
>
> In my case all ions are listed in the end of the pdb file, hence all IDs
> that are less than 344526 are protein atoms that I want to be left
> untouched. The force is applied to the rest of the ions if they enter the
> sphere of the radius and centre specified in the main configuration file.
>
> Hope that helps.
>
> Olga
>
>
> On Thu, Jan 12, 2017 at 3:42 AM, Azadeh Alavizargar <
> azadeh.alavizargar_at_gmail.com> wrote:
>
>>
>> Hello everyone
>>
>> I have a protein-membrane system with lots of ions. I want to prevent the
>> ions to enter a hypothetical cylinder aligned with the pore of the protein
>> using the upperBoundary of the distance xy of the colvar module. Do I have
>> to have a colvar for each ion or there is a way to do it for all the ions
>> with a single colvar?
>> If not, what is the best way to do it?
>>
>> Thank you very much for your help.
>>
>> Kind regards,
>> Azadeh
>>
>
>
This archive was generated by hypermail 2.1.6 : Sun Dec 31 2017 - 23:20:59 CST