vmd教程 forces-tutorial

vmd教程 forces-tutorial

UniversityofIllinoisatUrbana-Champaign

BeckmanInstituteforAdvancedScienceandTechnology

TheoreticalandComputationalBiophysicsGroup

ComputationalBiophysicsWorkshop

User-De nedForcesin

NAMD

AlekAksimentiev

DavidWells

GregSigalov

November2006Acurrentversionofthistutorialisavailableat

http://www.ks.uiuc.edu/Training/Tutorials/

Jointhetutorial-l@ks.uiuc.edumailinglistforadditionalhelp.

vmd教程 forces-tutorial

CONTENTS2Contents

1TclForces

1.1Introduction...............

1.2Example1:ConstantForces......

1.3Example2:Rotation.......................................................5559

1.4Example3:ForcingaSubsetofAtoms...............

1.5Example4:ImprovingE ciency..................

2TclBC

2.1Example1:MakingaBubble....................

2.2OptimizingaTclBCscript......................

2.3Example2:Concentratingtheionsinasolution.........

2.4Example3:Improvinge ciency..................

2.5Example4:ImposingaShearFlow.................

3Grid-steeredMolecularDynamics

3.1Introduction..............................

3.2Example1:ConstantForce.....................

3.3Example2:Surface..........................

3.4Example3:Multiplegrids......................13162021242731353939394346

vmd教程 forces-tutorial

CONTENTS3Introduction

ThistutorialisdesignedtoguideusersofVMDandNAMDintheuseoftclForcesandtclBC.Itisassumedthatyoualreadyhaveaworkingknowl-edgeofVMDandNAMD,aswellastheTcllanguage.FortheaccompanyingVMDandNAMDtutorialsgoto:

http://www.ks.uiuc.edu/Training/Tutorials/

Thistutorialhasbeendesignedspeci callyforVMD1.8.5,andshouldtakeabout5hourstocompleteinitsentirety.

Thetutorialissubdividedintotwoseparateunits,onefortclForcesandonefortclBC.Thesescript-basedfacilitiesmakepossiblevirtuallyanyperturbationofyoursystem,andaftercompletingthistutorialyouwillhavethebasicskillstousethem.

Throughoutthetext,somematerialwillbepresentedinseparate“boxes”.Theseboxesincludecomplementaryinformationtothetutorialandtipsortechnicaldetailsthatcanbefurtherexploredbytheadvanceduser.

Tclcodeiswritteninthisfont.Whenalineendswithabackslash(\)thismeansthatthereshouldbenolinebreak...don’tpressenter,justkeeptyping!Ifyouhaveanyquestionsorcommentsonthistutorial,pleaseemailtheTCBTutorialmailinglistattutorial-l@ks.uiuc.edu.Themailinglistisarchivedathttp://www.ks.uiuc.edu/Training/Tutorials/mailing

vmd教程 forces-tutorial

CONTENTS4Requiredprograms

Thefollowingprogramsarerequiredforthistutorial:

VMD:Availableathttp://www.ks.uiuc.edu/Research/vmd

NAMD:Availableathttp://www.ks.uiuc.edu/Research/namd

GettingStarted

Youcan ndthe lesforthistutorialintheforces-tutorial-filesdirectory.Fig.1showsthe lesanddirectoriesofforces-tutorial-files.

tclForcesFilescommonexample-outputpush

rot-a

rot-b

rot-c

BUBBLE

BUBBLE_OUTPUT_EXAMPLE

IONS

IONS_OUTPUT_EXAMPLE

SHEAR

SHEAR_OUTPUT_EXAMPLE

example-output

push-grid

graphene

ions-gridforces-tutorial-filestclBCfilestclBCmoviesgridForceFiles

Figure1:Directorystructureofforces-tutorial- les

TostartVMDtypevmdinaUnixterminalwindow,double-clickontheVMDapplicationiconintheApplicationsfolderinMacOSX,orclickontheStart→Programs→VMDmenuiteminWindows.

vmd教程 forces-tutorial

1TCLFORCES51TclForces

Inthisunit,youwilllearnaboutNAMD’sTclForcesfunctionality.AlthoughmanyuserneedswereanticipatedinthedevelopmentofNAMD,toanticipatethemallspeci callywouldbeimpossible.TclForcesallowsuserstoeasilysup-plementthebuilt-inconstraintandforcefunctionalitywithscript-basedforces.

1.1Introduction

ThebasicideaoftclForcesissimple:weprovideNAMDwithaTclscriptwhichtellsNAMDtoapplycertainforcestocertainatoms.Theonlyspecialrequirementofthisscriptisthatitde neacommandnamedcalcforces.ThiscommandiscalledbyNAMDeachtimestep.Withinthescript,wecanaccessatomicIDs,positionsandmasses.Inthisunit,wewillseehowtousetclForcesthroughfourexamplesofincreasingcomplexity.

1.2Example1:ConstantForces

Forour rstexposuretotclForces,wewilldosomethingverysimple:applyaconstantnetforcetoasystem.

1First,let’slookattheNAMDcon guration le.Gotothedirectoryforces-tutorial-files/tclForcesFiles/pushandopenthe lepush.namd.Everythingisstandardaboutituntiltheend,whereweseesomeunfa-miliarlines:

tclforceson

setlinaccel"3000"

tclforcesscriptpush.tcl

The rstlinesimplyturnstclForceson.Thesecondsetsaparameterfortheaccelerationthatwe’lluseinthescript,aswe’llsoonsee.Thelastlineisthenameofthescriptwhichactuallydoesthe

work.In-linescripts.Itisalsopossibletoincludethecontentsofyour

tclForcesscriptdirectlyinyourNAMDcon guration le.Instead

ofprovidingthenameofascript,wecouldinsteadwritethescript

contentswithincurlybraces:

tclforcesscript{

...contentsofscript...

}

Thisismostusefulforparticularlyshortscripts.

2Nowlet’slookatthescript.Openthe lepush.tcl,inthesamedirectory.

vmd教程 forces-tutorial

1TCLFORCES

3First,weselecttheatomswe’reinterestedin:

setnumatoms1231

setatoms{}

for{seti1}{$i<=$numatoms}{incri}{

lappendatoms$i

}

foreachatom$atoms{

addatom$atom

}6

Thestructurewe’reusingisubiquitin,thesameproteinusedintheVMDandNAMDTutorials.Here,we’reusingtheproteinalone,invacuum,andsothestructurecontainsjust1231atoms.We rstbuildalistoftheatomicindices(NAMDatomindices,unlikeVMDatomindices,startat

1.)Foreachatom,wethencalltheaddatomcommand.ThistellsNAMDthatwewantaccesstothisatom’scoordinates.Forcesmaybeappliedtoanatomwithoutthiscall.

4Wenextconvertfromthenaturalaccelerationunits A·ps 2totheNAMDunitskcal/mol· A·amu,andprintsomeinformation:

setlinaccel[vecscale[expr1.0/418.68]$linaccel]

print"Linearaccelerationimparted:($linaccel)Ang*ps^-2"Notethattoprintmessagesfromtheforcescript,wemustusetheprintcommand,insteadoftheusualputs

command.NAMDUnits.The

basicNAMDunitsarekcal/molforenergy,

angstromsforlength,atomicmassunits(a.k.a.Daltons)formass,

andbarforpressure.Otherunitsarederivedfromthese,e.g.the

NAMDunitofforceis1kcal/mol· A≈69.48pN.Linearacceleration.LinearaccelerationisgovernedbyNewton’s

SecondLaw:

F=ma

Here,Fisaforcevector,misthemass,andaistheacceleration

vector.Thus,thisequationsaysthattoachieveaspeci edaccel-

eration,wemustapplyaforceequaltothataccelerationtimesthe

mass,inthesamedirectionwewanttoaccelerate.Inourcase,we

applyaforcetoeachindividualatom,scaledbyitsmasssothat

eachexperiencesthesameacceleration,andthustheproteinmoves

togetherasawhole.

5Theonlythingleftisthecalcforcescommand:

vmd教程 forces-tutorial

1TCLFORCES

proccalcforces{}{

globalatomsnumatomslinaccelloadcoordscoords

loadmassesmasses

setcomsum"0.00.00.0"

settotalmass0.0

foreachatom$atoms{

setforce[vecscale$masses($atom)$linaccelnamd]

addforce$atom$force

settmp[vecscale$masses($atom)$coords($atom)]

setcomsum[vecadd$comsum$tmp]

}

setinvmass[expr1.0/$totalmass]

print"Centerofmass=[vecscale$invmass$comsum]"

}7

Thescriptstartswiththeglobalcommand,whichimportsthevariablesnamedsothattheycanbeaccessedinthecalcforcescommand.

Nextwecallthecommandloadcoords.ThissetsthevariablecoordstoaTclarrayofcoordinatesoftheatomsthathavebeenaddedusingaddatom.Inthiscase,ofcourse,thatisalltheatomsofthesystem.ArrayelementsinTclareaccessedusingparentheses.Similarly,loadmassesretrievestheatomicmasses,puttingtheminanarraynamedmasses.

Thenweloopthroughtheatoms,applyingthespeci edforcetoeachusingtheaddforcecommand.Alongtheway,wealsocalculatethecenterofmassoftheprotein.WeseetwomoreTclcommandsinactionshere:vecaddandvecscale.Theseaddtwovectorsandmultiplyavectorbyascalar,respectively.

6Becausethesystemissosmall,youcanrunitonyourdesktoporlaptopwithouttrouble,andseetheproteinmove.Runthe lepush.namd.Withanormalinstallation,youwoulddothiswiththecommand

namd2push.namd>

push.logCoordinates.CoordinatesintclForcesarenotwrappedaround

theunitcell.Ifyouneedwrappedcoordinatesforwhateverreason,

youmustwrapthemyourselfusingtheperiodiccelldimensionsand

origin.

7Afterthis nishes,wewanttoseehowtheproteinwasa ected.WewillnowuseVMD’sMultiPlotplugintoviewtheresults.Openthe leplot.tclinyourfavoritetexteditor,andchangethe rstlinetobethelocationofthelog leyoujustcreated,e.g.

setlogpush.log

ifyounamedthelog lepush.log.Torunthescript,eithertype

vmd教程 forces-tutorial

1TCLFORCES

vmd-eplot.tcl8

fromthecommandline,oropenVMDandsourcethescript leintheTkConsole:

sourceplot.tclCenter of mass

70

60

50

x (Å)40

30

20

10

00500100015002000time (ns)25003000

Figure2:Centerofmassposition

plottedversustime.

8YoushouldseeaplotsimilartoFig.2.Thebeginningofthegraphisquadratic,asweexpect.However,ngevindy-namicsinvolveadampingfactor,ngevindynamicsarecommonlyusedin

NAMDtosimulateatconstanttemperature,inthiscaseinaso-

calledNVTensemble.Thekeywordsbeginningwithlangevinin

theNAMDcon guration lecontrolitssettings.WithoutLangevin

dynamicsactivatedandintheabsenceofanyexternalforces,we

wouldbesimulatingintheNVEensemble,inwhichcaseenergy

ratherthantemperatureisconstant.NAMDalsoincludesconstant-

pressureoptionsforsimulationintheNPTensemble.

9NowopenthetrajectoryinVMD.FirstlaunchVMD,thenopentheTkConsole.Makesureyou’reinthetclForcesFilesdirectory,thentypethecommand

molloadpsfcommon/ubiquitin.psfdcdpush/push.dcd

Ifyoucouldnotproducethetrajectory,thereisasampleoneintheexample-outputdirectory.NoticethattheLangevindynamicspreventtheproteinfrommovingtooquickly.

vmd教程 forces-tutorial

1TCLFORCES9

1.3Example2:Rotation

Wewillnowaddalayerofcomplexity:wewillforceatomsdi erentlybasedontheircoordinates.Inthisexample,wewillforceatomsinsuchawaythat,inadditiontoalinearforce,theproteinalsorotatesaboutanaxisparalleltothez-axisandthroughtheprotein’scenterofmass.

1Changetothedirectoryforces-tutorial-files/tclForcesFiles/rot-a.2LookattheNAMDcon guration le.Openthe lerot-a.namd.Attheend,weseethefollowing:

tclforces

setlinaccel

setangaccel

tclforcesscripton"3000"1rot-a.tcl

Thisisverysimilartothepreviousexample,buthereweprovideanadditionalscalartodescribetheangularacceleration.

3Nowlet’slookatthescriptitself.Openthe lerot-a.tcl.Theatomselectionpartatthebeginningisidenticaltothelastcase,sowewillnotrepeatthediscussion.Thenextpartofthe leisagainsimilar,butalsoprocessestheangularacceleration.

setlinaccel[vecscale[expr1.0/418.68]$linaccel]

setangaccel[expr$angaccel/418.68]

print"Linearacceleration:($linaccel)Ang*ps^-2"

print"Angularacceleration:(00$angaccel)Rad*ps^-2"

WeconverttoNAMDunits,justlikebefore.Thenweprintsomeinfor-mationaltextusingtheprintcommand.

4Nowwe’lllookatthecalcforcescommand.The rstpartisnearlythesameasinthepreviousexample.Thedi erencesarethechangeofoneoftheglobalvariables,andthefactthatwenolongerapplyforcesinthe rstloop.

proccalcforces{}{

globalatomsnumatomslinaccelangaccelnamd

loadcoordscoords

loadmassesmasses

vmd教程 forces-tutorial

1TCLFORCES

setcomsum"000"

settotalmass0

foreachatom$atoms{

settmp[vecscale$masses($atom)$coords($atom)]

setcomsum[vecadd$comsum$tmp

settotalmass[expr$totalmass+$masses($atom)]

}

setcom[vecscale[expr1.0/$totalmass]$comsum]

print"Centerofmass=$com"

5Thelastpartofthescriptiswheretherealdi erenceslie:

foreachatom$atoms{

setlinforce[vecscale$masses($atom)$linaccelsetr[vecsub$coords($atom)$com]

setx[lindex$r0]

sety[lindex$r1]

setrho[exprsqrt($x*$x+$y*$y)]

setphi[expratan2($y,$x)+$MPI/2]

if{$atom==1}{

print"atom$atom:

}phi=$phi"10

setangdir"[exprcos($phi)][exprsin($phi)]0.0"

setangmag[expr$masses($atom)*$angaccelnamd*$rho]

setangforce[vecscale$angmag$angdir]

setforce[vecadd$linforce$angforce]

addforce$atom$force

}

}

The rstsectionabovecalculatesthedistancefromtherotationaxis,rho,andtheanglefromthexaxisatwhichtheforceshouldbeapplied,phi(90 fromtheangleofthepositionvector,hencetheadditionofπ/2radians.)Fig.3showshowthesevariablesarerelated.

Afterprintingtheforceangleassociatedwithoneatomsowecanmonitortherotation,we ndaunitvector(i.e.avectorwithlengthone)pointinginthephidirection.Wethenmultiplythisbythemassoftheatomanditsdistancefromourrotationaxis,givingusthecorrectangularforcefortheangularaccelerationrequestedintheNAMDcon guration le.Finally,weaddthelinearandangularforcestogetthetotalforce.

vmd教程 forces-tutorial

1TCLFORCES11x

Figure3:Relationshipsamongvariablesinthescriptrot-a.tclAngularacceleration.Angularaccelerationisgovernedbyan

equationanalogoustoNewton’sSecondLaw:

τ=Iω(1)

τistorque(analogoustoforce),ωistheangularacceleration,and

Iisthemomentofinertia(analogoustomass).Torqueisgivenby

τ=ρ×F

whereFistheforceandρisapositionvectorrelativetotherotation

axis(seeFig.3).Inourcase,weapplytheforceperpendicularto

theρvector,sothecrossproductsimpli esandweget

τ=ρF(2)

forthemagnitude.Forpointparticles,themomentofinertiais

simplyI=mρ2,wheremistheparticle’biningthis

withEqns.1and2andsolvingforF,weget

F=mωρ

fortheforceperparticletoachieveaspeci edangularacceleration.

vmd教程 forces-tutorial

1TCLFORCES

6Nowrunthesimulation:

namd2rot-a.namd>rot-a.log

Thisshouldagaintakejustafewminutes.12

7Plottheanglephibyusingthescriptplot.tcljustasyoudidbefore(thereisadi erentcopyinthecurrentdirectory.)

TheresultisshowninFig.4.Plotthecenterofmassaswell,byopeningplot.tclandchangingthevarvariableto"Centerofmass".Thiscanbeaccomplishedbysimplycommentingthe rstsetvarline,anduncommentingtheother:

#setvar"atom1:phi"

setvar"Centerofmass"Angle of rotation

2

1.5

1

0.5

-0.5

0φ (rad)50010001500time (fs)200025003000

Figure4:Rotationangleplottedversustime.

8Weagainseearoughlyquadratictimedependenceatthebeginningofthetrajectorygivewaytoalinearone,againaconsequenceoftheLangevindamping.

9OpenthetrajectoryinVMD.WithVMDopen,gotothemaintclForcesFilesdirectory,thentypethefollowingintheTkConsole:

molloadpsfcommon/ubiquitin.psfdcdrot-a/rot-a.dcd

Challenge:Findthedependenceoftheterminallinearandangularveloci-tiesonthelangevinDampingparameterofthesimulation,setinthecon gura-tion le,aswellasthedependenceontheappliedforce.

vmd教程 forces-tutorial

1TCLFORCES13

1.4Example3:ForcingaSubsetofAtoms

Thisexampledealswiththeimportantabilitytoforceasubsetofyouratoms.Inthissection,wewillapplyforcesonlytothebackboneatomsoftheprotein.Weaccomplishthisbymakingaspecial“targetatom”PDB le,usingthebetacolumntomarktheatomswewanttoapplyforceto.Inaddition,wewillusetheoccupancycolumntotellourscripttheatomicmasses,inordertodemon-stratehowtoprovidethescriptwithadditionalinformation.

1Changetothedirectoryforces-tutorial-files/tclForcesFiles/rot-b.2Asusual,wewill rstexaminetheNAMDcon guration le.Openthe lerot-b-short.namd.First,notethatwearenowusingadi erentsystem,oneinwhichtheubiquitinhasbeensolvated.Thisisamorenaturalenvironmentfortheprotein.

Lookattheendofthe le,wherewesetuptclForces.Weseeanewline:

settargetAtomPdb../common/ubiquitinbackbone.pdb

Withthisline,wespecifytheaforementionedtargetatomPDB le.Thus,our rsttaskistoproducethat le.WewilluseaVMDscripttoaccom-plishthis.

3Openthe lemakeTargetAtomPdb.tcl.The rstsectionsetsafewvari-ables.First,therelevant lenames:

setpdbcommon/ubiquitinsetpsfcommon/ubiquitinsettargetPdbcommon/ubiquitinNext,wesettheselectiontextwewillusetoselectouratoms,andthebetavaluethattheseatomswillreceive:

setselection"proteinandbackbone"

settargetMark"1.00"

4Next,weloadthestructureintoVMD,andsetboththebetaandoccu-pancycolumnstozero:

molloadpsf$psfpdb$pdb

setall[atomselecttopall]

$allsetbeta0

$allsetoccupancy0

5Nowwemakeaselectionofthetargetatoms,settheirbetaandoccupancycolumnstotheappropriatevalues,andwritetheoutputPDB.Notetheeasewithwhichthelistofmassesareusedtosettheoccupancy:

vmd教程 forces-tutorial

1TCLFORCES

settarget[atomselecttop$selection]

setmasses[$targetgetmass]

$targetsetbeta$targetMark

$targetsetoccupancy$masses

$allwritepdb$targetPdb

exit

6Nowrunthescript:

vmd-dispdevtext-emakeTargetAtomPdb.tcl14

(Note:ForWindowsusers,opentheVMDGUIandsourcethescript leintheTkConsole.)

Ifyouwereunabletorunthescript,copythe leexample-output/ubiquitinsolvateintothecommondirectory.

7Next,takealookatthetclForcesscript.Openthe lerot-b.tcl.Rightfromthebeginning,wenoticemanydi erences.It rstsetsatargetMarkvariablesothatthescriptcanrecognizethemarkswe’vesetinthetargetPDB:

settargetMark1.0

8NowwehavetoprocessthePDB.Therearenobuilt-inroutinesforac-complishingthis,sothefollowingcodeisquitelow-level.ThecolumnsofaPDB lehavea xedcharacterwidth,sowemustsimplyreadeachline,andbreakitinto xed-sizedpiecesaccordingtothePDBformat:

settargets{}

setmasses{}

setinStream[open$targetAtomPdbr]

foreachline[split[read$inStream]\n]{

settype[stringtrim[stringrange$line05]]

setname[stringtrim[stringrange$line1215]]

setresid[stringtrim[stringrange$line2225]]

setbeta[stringtrim[stringrange$line6065]]

setoccupancy[stringtrim[stringrange$line5459]]

setsegname[stringtrim[stringrange$line7275]]

9We rstmakesurethatthislinecorrespondstoanatomrecordandnotacommentlineorsomeotherentry.Then,iftheatomhasabetavaluematchingthetargetvalue,weformatripleconsistingoftheatom’ssegname,resid,andname.Thisisnecessaryfor ndingtheindexoftheatom:

vmd教程 forces-tutorial

1TCLFORCES

if{($typeeq"ATOM"||$typeeq"HETATM")\

&&$beta==$targetMark}{

lappendtargets"$segname$resid$name"

lappendmasses$occupancy

}

}

close$inStream15

10Thenextstepistousetheatomtriplesto nditsindex,andformalist

ofthese,usingtheatomindexcommand.addatomisthencalledforeachoftheseatoms:

setatoms{}

foreachtarget$targets{

lassign$targetsegnameresidatom

setatomindex[atomid$segname$resid$atom]

lappendatoms$atomindex

addatom$atomindex

}

11Nextwe ndthenumberoftargetatoms:

setnumatoms[llength$atoms]

setlinaccel[vecscale[expr1.0/418.68]$linaccel]

setangaccel[expr$angaccel/418.68]

print"Linearacceleration:($linaccel)Ang*ps^-2"

print"Angularacceleration:(00$angaccel)Rad*ps^-2"

12Nowwegettothemainpartofthescript,thecalcforcesde nition.It

isessentiallyidenticaltothelastexample:

proccalcforces{}{

globalatomsnumatomsmasseslinaccelangaccelnamd

loadcoordscoords

setcomsum"000"

settotalmass0

foreachatom$atomsmass$masses{

settmp[vecscale$mass$coords($atom)]

setcomsum[vecadd$comsum$tmp]

settotalmass[expr$totalmass+$mass]

}

setcom[vecscale[expr1.0/$totalmass]$coordsum]

print"Center=$com"

vmd教程 forces-tutorial

1TCLFORCES

foreachatom$atomsmass$masses{

setlinforce[vecscale$mass$linaccelnamd]

setr[vecsub$coords($atom)$com]

setx[lindex$r0]

sety[lindex$r1]

setrho[exprsqrt($x*$x+$y*$y)]

setphi[expratan2($y,$x)+$MPI/2]

if{$atom==1}{

print"atom$atom:

}phi=$phi"16

setangdir"[exprcos($phi)][exprsin($phi)]0.0"

setangmag[expr$angaccel*$rho*$mass]

setangforce[vecscale$angmag$angdir]

setforce[vecadd$linforce$angforce]

addforce$atom$force

}

}

13Nowrunthesimulation,justasyouhaveinpreviousexamples:

namd2rot-b-short.namd>rot-b-short.log

Ifyouhaveaccesstoacluster,runitonthecluster—refertolocalin-structionsforhowthisisdone.Thissimulationwillrunforjust100steps,sothatthebenchmarktimesareprinted.Inthenextsection,wewillseehowtoimprovethisperformance.

Ifyouareinterested,examinethe lesrot-b-long.logandrot-b-long.dcdintheexample-outputdirectory,whichwereproducedwithaslightlylongersimulation.WhenloadingthetrajectoryinVMD,usethePSFcommon/ubiquitinsolvate.psf

1.5Example4:ImprovingE ciency

Inmostcases,it’sunnecessarytorecalculatetheforcevaluesappliedateachtimestep.Dependingonthesituation,onemaybeabletorecalculateevery1000steps,savingalotofcomputationale ortandthereforeincreasingthespeedofthesimulation.Inthissection,wewillseehowtodothat.

1Changetothedirectoryforces-tutorial-files/tclForcesFiles/rot-c.2OpentheNAMDcon guration le,rot-c-short.namd.Thenewlinethistimeis

setforcesRecalcFreq10

Asitsnameimplies,thisparametersetshowoftenourforceswillberecalculated.

vmd教程 forces-tutorial

1TCLFORCES173Nowopenthescriptrot-c.tcl.The rstdi erencesstartatline58.We rstsetupaforceslist,whichisavariabletoholdthevaluesoftheforceswewillapplytotheatoms.Wealsosetupsomecountervariables:forcecountandprintcountarebothincrementedeachtimestep.WhenforcecountequalsforcesRecalcFreq,theforcesarerecalculated,andforcecountisresettozero,asseenbelow.

setforces{}

foreachindex$atoms{

lappendforces"0.00.00.0"

}

setforcecount$forcesRecalcFreq

setprintcount0

4Nowwebeginthecalcforcescommand.First,asusual,wedeclareourglobalvariables.Wethenapplytheforcesthataresavedintheforceslist.

proccalcforces{}{

globalatomsnumatomsforcemultmassesavgmassforces

globalforcesRecalcFreq

globalforcecountprintcount

foreachatom$atomsforce$forces{

addforce$atom$force

}

5Next,wetestwhethertheforceswillberecalculatednexttimestep,andifso,tellNAMDthatwewillwantatomiccoordinatesforthetargetatomsbycallingaddatom:

if{$forcecount==[expr$forcesRecalcFreq-1]}{

print"Addingatomspriortoreconfiguringforcesat\

$printcount"

foreachatom$atoms{

addatom$atom

}

}

Aswewillseeshortly,wecallclearconfigaftertheforcesarerecalcu-lated.Thiscallerasesalladdatomrecords;withoutit,thecoordinatesoftheatomsaddedwillbeavailableeverytimestep,independentlyofwhetherwecallloadcoords,andthereforemuchofthepotentialspeedgainwillbelost.However,thismeansthataddatommustbecalledagaineachtimewewanttorecalculateforces,andbecauseoftechnicaldetailsithastobedoneatleastonestepbeforethecoordinateswillbeused.

6Nextwehavethecodethatrecalculatestheforce.Asalludedtoabove,therecalculationhappenswhenforcecountequalsforcesRecalcFreq.

vmd教程 forces-tutorial

1TCLFORCES18

Therestoftheforcecalculationcodeisidenticaltothecodewehadinthelastexample,exceptthatinsteadofcallingaddforcewitheachforcevector,wenowputthevectorintheforceslistvariable:

if{$forcecount==$forcesRecalcFreq}{

print"Recalculatingforcesat$printcount"

loadcoordscoords

setcomsum"000"

settotalmass0

foreachatom$atomsmass$masses{

settmp[vecscale$mass$coords($atom)]

setcomsum[vecadd$comsum$tmp]

settotalmass[expr$totalmass+$mass]

}

setcom[vecscale[expr1.0/$totalmass]$comsum]

print"Centerofmass=$com"

setforces{}

foreachatom$atomsmass$masses{

setlinforce[vecscale$mass$linaccelsetr[vecsub$coords($atom)$com]

setx[lindex$r0]

sety[lindex$r1]

setrho[exprsqrt($x*$x+$y*$y)]

setphi[expratan2($y,$x)+$Mif{$atom==1}{

print"atom$atom:

}phi=$phi"

setangdir"[exprcos($phi)][exprsin($phi)]0.0"

setangmag[expr$angaccel*$rho*$mass]

setangforce[vecscale$angdir]

setforce[vecadd$linforce$angforce]

lappendforces$force

}

7Finally,weprintsomeinformation,callclearconfig,andresetforcecount.

print"Step${printcount}:Recalculated\

[llength$forces]forces"

setforcecount0

clearconfig

}

incrforcecount

incrprintcount

return

}

vmd教程 forces-tutorial

1TCLFORCES

8Nowrunthesimulation:

namd2rot-c-short.namd>rot-c-short.log19

Again,thisisashortsimulation,running100steps,parethespeedofthissimulationtothelast.Thereshouldbeasigni cantdi erence—comparingtheexample lesrot-b-short.logandrot-c-short.log,whichwererunonjustasingleCPU,thereisnearlya15%speedincrease.AsthesizeofthesystemandthenumberofCPUsincrease,applyingtclForcesef- cientlybecomesevenmoreimportant:comparingrot-c-long.logwithrot-b-long.log,bothrunon20processors,rot-cisalmost60%faster!

Onceagain,therearesampleoutput les,rot-c-long.dcdandrot-c-long.logintheexample-outputdirectory,andshouldagainbeloadedtogetherwiththePSF lecommon/ubiquitinsolvate.psf.

Challenge:Taketheprovided lekcsa,andusethepartialstructurekcsaasatargettoopenthechannelusingtclForces.YounowknowthebasicsoftclForces.Inthenextsection,youwilllearnaboutasimilarfeatureofNAMDcalledtclBC,whichissuitablefordi erentcircumstances.

vmd教程 forces-tutorial

2TCLBC202TclBC

InaNAMDcon guration le,onesetsupthesimulationparameters,suchasthetemperature,pressure,orauniformelectric eld,forthemolecularsystemasawhole.Inmanycases,though,youwillneedtoimposeboundaryconditions,applyanon-uniformelectric eld,orselectivelyapplyforcestoagroupofatomsofthegiventypeorlocatedinaparticulararea.Forthispurpose,youcanuseatclBC(whichstandsforTclBoundaryConditions)script.AtclBCscriptiswrittenintheTclprogramminglanguage,arudimentaryknowledgesofwhichwouldbeenoughformosttasks.Thescriptwillbeinterpretedline-by-line,thereforeslowingdownthesimulation,sometimesconsiderably.Forthisreasonyouwillwanttokeepyourscriptssimpleand

e cient.tclBCversustclForces.tclBCandtclForcesarebothscripting

interfacesthatallowtogetinformationabouttheatoms’positions

andapplyforcestoselectedatomsduringaNAMDsimulation.In

bothscripts,theforcecanbecalculatedindividuallyforeachatom,

dependingonitsuniqueIDnumberanditscurrentposition.The

maindi erencebetweentclBCandtclForcesisthatanindepen-

dentinstanceoftclBCisrunningoneachprocessor,andonlyatoms

formingthepatchtreatedbythatprocessorarevisibletothegiven

instanceoftclBC.ThisfeaturemakestclBCmoree cientthan

tclForces,whichisrunonjustoneCPU,butalsolimitsitsca-

pabilities.Ifyouneedtoapplyforcestoeachatomirrespectiveof

thepositionofotheratoms,usetclBC.Ontheotherhand,ifyou

havetoconsidermutualpositionsoftwoormoreatomsorgather

informationaboutthewholesystem,youwillneedtclForces.

ThewaythecoordinatesofatomsarewrappedaroundperiodicboundariesinatclBCscriptiscontrolledusingthecommandwrapmode,whichthatmayhaveoneofthefollowingarguments:

patchisthedefaultmode:theatom’scoordinatesareconsideredasthepositioninNAMD’sinternalpatchdatastructure.Thisisingeneraldi erentfromanatom’scoordinatesrelativetotheglobalorigin,sopatchmodeshouldnotbeusedunlessyouknowwhatyouaredoing.

Inmodeinput,theatom’scoordinatescorrespondtoitspositionintheinput lesofthesimulation.

Inmodecell,theatom’sequivalentpositionintheunitcellcenteredonthecellOriginisused.

Inmodenearest,theatom’sequivalentpositionnearesttothecellOriginisused.

Formostpurposes,modecellisagoodchoice.

vmd教程 forces-tutorial

2TCLBC21

2.1Example1:MakingaBubble

Imaginethatyouhaveawaterbox,andyouwanttocreateasphericalbubbleofvacuum.Youcandothatbyapplyingforcestoeveryatomfoundinsidethisspheretopushitout.Toavoidinstability,youwouldwanttostartwithasmallbubbleandthenincreaseituntilthedesiredsizeisreached.Therateofincreasingthesizeofthebubbleisuptoyou.Ifthisrateistoolarge,someatomswillbepushedtoohardandthereforemovetoofast,causingthesimulationtocrush.Fast-movingatomsmayalsobreakthestructureofadouble-strandedDNAoranothermolecularcomplexheldtogetherbyrelativelyweaknon-bondedinteractions(vanderWaalsandelectrostaticforces).Evenifthesimulationremainsstable,youshouldwatchthechangesintemperatureandpressuretomakesurethattheenergyin ux,whichisequaltotheworkdonebytheforceyouareapplying,hasenoughtimeto

dissipate.Choosingtheloadingrate.Generallyyouwouldwanttokeep

theperturbationofthesystemaslowaspossiblebyimposingany

externalforcesasslowlyasyoucana ordgivenyouravailablecom-

putationalresources.Itisagoodidea,though,torunatrialsim-

ulationatafastpacetoseewhatishappeningtothesystem,and

therebygetanideaofhowfastyoucanmovethingswithoutlos-

ingthestabilityofthesimulationorincreasingthetemperatureor

pressurebeyondreason.

1Openthe letclBCfiles/BUBBLE/eq04.bubble.namdinatexteditor.Scrolldowntotheendofthe letothelinetclBCon.

Youshould ndthefollowingcode:

tclBCon

tclBCScript{

setbubbleCenter"0.00.00.0"

settclBCScriptbubble.tcl

source$tclBCScript

}

tclBCArgs{0.15.0.015.}

SinceatclBCscriptiscalledfromNAMD,itisreferencedNAMDcon- guration le,whichisalsoagoodplacetosetthescriptparameters.TheNAMDcommandtclBConturnstheTclBCinterfaceon.Inourexample,tclBCScript{...}containstheinitializationofakeyvariableandareferencetothe lethatcontainsthescriptitself:source$tclBCScript.Ifthescriptisshort,itcanbeplacedentirelywithinthebodyofthecommandtclBCScriptintheNAMDcon guration le.Finally,thecommandtclBCArgsisusedtopassalistofvariablestothemainTclBCprocedurecalcforces(wewilltalkaboutitverysoon).Inthiscase,thefourargumentsfoundincurlybracketshavethemeaning:“makeabubblestartingfromradius0 Aandincreaseitto15 Aata rateof0.01Apersimulationstep,byapplyingforcesof5kcal/mol·A.”

本文来源：https://www.bwwdw.com/article/qp3j.html