From: Hongbo Du (hongbodu_at_uark.edu)
Date: Fri Feb 28 2014 - 11:22:04 CST
Hi, NAMD developers and users,
I am trying to simulate salt solutions with the external magnetic field through TCLBC or TCLforces. I have read out the previous and current coordinates of all the atoms and I want to use the equation F=q*(V cross B) to add a magnetic force on each atom. I specify the magnetic intensity of 0.05 T or 10 T and found that there almost doesn't exist effects on the salt solutions. I tried to use the commands loadtotalforces and loadforces available in TCLforces to check the forces I applied and the internal forces. It was found that the applied force was tiny and could be neglected. The reason is that the number of the velocity is small when velocity is obtained from the differences between the current and previous coordinates divided by timestep. Another way, Can I read out the velocity through TCLBC directly? Do you think there are some problems I may have when I add the magnetic force through TCLBC. The following is my code through TCLBC.
##################
wrapmode cell
set numatoms 6710
# Two first agruments of calcforces are automatically forwarded
# to it by NAMD. The other 4 arguments match the list of 4 values
# from command tclBCArgs.
for { set i 1 } { $i <= $numatoms } { incr i } {
set previouscoorx($i) 0
set previouscoory($i) 0
set previouscoorz($i) 0
}
#for { set i 1 } { $i <= $numatoms } { incr i } { print "previous $previouscoor($i)" }
proc calcforces { step unique } {
global forceunit charge mfield diffcoord atomnumber
global previouscoorx previouscoory previouscoorz
set frametime 0.000000000000002
set forceunit 6.948E-11
set mfield "0.05 0 0"
if { $step == 1 } {
while {[nextatom]} {
set atomnumber [getid]
set rvec [getcoord]
foreach {x y z} $rvec { break }
set previouscoorx($atomnumber) $x
set previouscoory($atomnumber) $y
set previouscoorz($atomnumber) $z
# print "$atomnumber test x $previouscoorx($atomnumber)"
}
# print "step $step"
}
if { $step >= 2 } {
while {[nextatom]} {
set atomnumber [getid]
set charge [getcharge]
set rvec [getcoord]
set tmpcharge [expr $charge*1.602176565E-19]
set mfieldx [lindex $mfield 0]
set mfieldy [lindex $mfield 1]
set mfieldz [lindex $mfield 2]
foreach {currentx currenty currentz} $rvec { break }
set currentvx [expr ($currentx-$previouscoorx($atomnumber))*1.0E-10/$frametime]
set currentvy [expr ($currenty-$previouscoory($atomnumber))*1.0E-10/$frametime]
set currentvz [expr ($currentz-$previouscoorz($atomnumber))*1.0E-10/$frametime]
set mforcex [expr $tmpcharge*[expr $mfieldy*$currentvz-$mfieldz*$currentvy]/$forceunit]
set mforcey [expr $tmpcharge*[expr -$mfieldx*$currentvz+$mfieldz*$currentvx]/$forceunit]
set mforcez [expr $tmpcharge*[expr $mfieldx*$currentvy-$mfieldy*$currentvx]/$forceunit]
addforce "$mforcex $mforcey $mforcez"
if { ($atomnumber==6710) && ( $step%100==0) } {
print "atom $atomnumber previous $previouscoorx($atomnumber) $previouscoory($atomnumber) $previouscoorz($atomnumber)"
print "atom $atomnumber current $currentx $currenty $currentz "
print "velocity $currentvx $currentvy $currentvz"
print "step $step atom $atomnumber $mforcex $mforcey $mforcez"
}
set previouscoorx($atomnumber) $currentx
set previouscoory($atomnumber) $currenty
set previouscoorz($atomnumber) $currentz
}
#print "step $step"
}
while {[nextatom]} {
}
}
########################
Thank you very much for your help in advance.
Best wishes,
Dr. Hongbo Du
Postdoc fellow at University of Arkansas
This archive was generated by hypermail 2.1.6 : Wed Dec 31 2014 - 23:22:11 CST