Re: Use tclforces to find systematic forces on fixed atoms?

From: Bill Kowallis (bill_at_mercury.chem.pitt.edu)
Date: Wed Oct 08 2008 - 11:45:12 CDT

Hi Jerome,
    Yes sir, I have tried the fixedAtomsForces set to off and on
with the same result. I'll paste my namd script above my tcl script
if you wish to take a look. Thank you very much for your response,
I appreciate it.
    -Bill

> Hi Bill,
> Are you using the fixedAtomsForces option?
> Jerome
>
> On Tue, Oct 7, 2008 at 5:20 PM, Bill K <bill_at_mercury.chem.pitt.edu> wrote:
>> Hello everyone,
>> I've encountered an interesting (and yet undiscussed?) problem using
>> the tclforces interface. I'm trying to keep a couple of atoms fixed
>> (in different areas of my simulation box - they're not interacting
>> with each other) and measure the forces exerted upon them by the
>> mobile atoms of the system. I used a fixedatomfile and set my target
>> atoms' beta values to 1 in the fixed atom pdb to immobilize them. My
>> tcl script works fine, by the way, as I've tested it on mobile atoms.
>> The problem is that calcforces proc prints out values of 0.0 0.0 0.0
>> for the force on these fixed atoms.
>> When I do this calculation, the fixed atoms must be completely
>> restrained, instead of simply applying a harmonic restraint, since
>> I'm calculating a force autocorrelation function. If anyone has
>> tackled this problem but just didn't post it, or has a good
>> suggestion, I would greatly appreciate your time and consideration.
>> I'm thinking that calcforces is printing out the forces after the
>> fixed atom forces have already been set to zero (which they must be
>> in order for their displacement to be zero). I'll paste my tcl
>> script below just in case you're curious. Thank you very much for
>> your time.
>> -Bill
>>

structure protNOV5-3.psf
coordinates free_prot_200mV_01.coor
velocities free_prot_200mV_01.vel
extendedSystem free_prot_200mV_01.xsc
parameters top_all27_prot_lipid_dum.inp
paraTypeCharmm on
outputname facf_prot_200mV_01

outputEnergies 200
outputTiming 200
dcdFreq 2000
restartFreq 2000
wrapAll on
wrapNearest on

timestep 3
rigidBonds all
nonBondedFreq 1
fullElectFrequency 2
stepsPerCycle 10

switching on
switchDist 7.5
cutoff 9
pairlistdist 10.5

margin 5

Pme on
PmeGridsizeX 96
PmeGridsizeY 96
PmeGridsizeZ 128

exclude scaled1-4
1-4scaling 1.0

fixedAtoms on
fixedAtomsForces on
fixedAtomsFile freeprot_freeze2ions.pdb
fixedAtomsCol B

langevin on
langevinDamping 1
langevinPistonTarget 1.01325
langevinPistonPeriod 200
langevinPistonDecay 100
langevinPistonTemp 310
useGroupPressure yes
useFlexibleCell no
useConstantRatio no
useConstantArea no

eFieldOn no
eField 0.0 0.0 0.036

binaryoutput off

tclForces on
tclForcesScript force_calculation.tcl

run 2000

>>
>>
>> set atm1 250
>> set atm2 4560
>>
>> addatom $atm1
>> addatam $atm2
>>
>> set Tclfreq 10
>> set t 0
>>
>> set outfilename forces1_tcl1.out
>> open $outfilename w
>>
>> proc calcforces {} {
>>
>> global atm1 atm2 Tclfreq
>> global t outfilename
>>
>> loadtotalforces forc
>>
>> if {$t>1} {
>> if {[array exists forc]} {
>> print "Forces success"
>> set f1 $forc($atm1)
>> set f2 $forc($atm2)
>>
>> print "$f1 $f2"
>>
>> set outfile [open $outfilename a]
>> set time [expr $t/1]
>> puts $outfile "$time $f1 $f2"
>> close $outfile
>> }
>> } else {
>> print "Forces Failure (after waiting for 1 timestep)"
>> }
>> incr t
>> return
>> }
>>
>>
>>
>
>

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:49:56 CST