一维分子动力学模拟python代码

1#!/usr/bin/envpython323#Moleculardynamicsin1D45#importneededmodules6importrandom7importmath8importsys910#themainloopfunction11defmain(md):12time=0.0#initialisetime1314#opencoordinateoutputfiles15cfile=open(coords.xyz,w)16#opentemperatureoutputfile17tfile=open(temperature.dat,w)18tfile.write(#timetemperature\n)19#openenergyoutputfile20efile=open(energy.dat,w)21efile.write(#timeenergy\n)2223#mainMDloop24fortinrange(md.tsteps):25print(#----Time=,round(time,2),----Steps=,t,----)#python326en=md.force()#calculateforces27md.integrate(t,en)#integrateequationsofmotion28#printcurrentcoordinatestofile29md.printcoords(time,cfile)30#printcurrenttemperaturetofile31tfile.write(str(round(time,2))++str(md.temp)+\n)32#printcurrentenergytofile33efile.write(str(round(time,2))++str(md.etot)+\n)3435time+=md.dt#increasetimebydt3637md.statistics(tfile,efile)#calculateaverages,SD,etc3839#closeoutputfiles40cfile.close()41tfile.close()42efile.close()4344classMD(object):4546N=36#numberofparticles(integerforloopcontrol)47dN=36.0#numberofparticles(doublefordoingmaths)48L=36.0#lengthof1Dbox49dim=1.0#dimensions50#initialisepositionsandvelocites51def__init__(self):5253#declaretheglobalvariables54#constants55self.a=self.L/self.dN#latticespacing56self.dtsteps=100.0#numberoftimesteps(doubleformath)57self.tsteps=int(self.dtsteps)#numberoftimesteps(integerforcounting)58self.dt=0.01#integrationtimestep59self.rc=18.0#distancecutoffforcomputingLJinteractions60self.rc2=self.rc**2#distancecutoffsquared61self.ecut=4.0*((self.rc**-12)-(self.rc**-6))#valueofLJpotentialatr=rc:4(1/rc^{12}-1/rc^{6})6263#valuesthatchangeduringthesimulation64self.en=0.0#potentialenergy65self.etot=0.0#totalenergy(pot+kin)66self.temp=0.728#temperature6768#listsforstoringstuff69self.x=[]#coordinates70self.xp=[]#previouscoordinates71self.v=[]#velocities72self.f=[]#forces7374#storedataforaveraging75self.sumTemp=0.076self.sumEtot=0.077self.temps=[]78self.etots=[]7980#initialisetheliststobethecorrectsize81foriinrange(self.N):82self.x.append(0.0)83self.xp.append(0.0)84self.v.append(0.0)85self.f.append(0.0)8687#initialisetheliststobeofthecorrectsizefordatastorage88foriinrange(self.tsteps):89self.temps.append(0.0)90self.etots.append(0.0)91print(#----Initialisingpositionsandvelocities----)9293sumv=0.0#sumofvelocities94sumv2=0.0#sumofvelocitiessquared9596#loopovertheparticles97foriinrange(self.N):98self.x[i]=self.lattice_pos(i)#placeparticlesonalattice99self.v[i]=random.random()-0.5#assignvelocitiesuniformly(butrandomly)inrange[-0.5,0.5]100sumv+=self.v[i]#sumvelocities101sumv2+=self.v[i]**2#sumsquardvelocities102103sumv=sumv/self.dN#velocityofcentreofmass104sumv2=sumv2/self.dN#mean-squaredvelocity105sf=math.sqrt(self.dim*self.temp/sumv2)#scalefactorforvelocitestoachievedesiredtemperature106107foriinrange(0,self.N):108self.v[i]=(self.v[i]-sumv)*sf#scalevelocites109self.xp[i]=self.x[i]-self.v[i]*self.dt#setpreviouspositions110111112#placeparticlesonalattice113deflattice_pos(self,i):114pos=(i+0.5)*self.a115returnpos116117#calculateforces118defforce(self):119print(#----Calculatingforces----)120en=0.0#(re)setenergytozero121#(re)setforcestozero122foriinrange(self.N):123self.f[i]=0.0124125#loop(inefficiently)overallpairsofatoms126foriinrange(0,self.N-1):127forjinrange(i+1,self.N):128xr=self.x[i]-self.x[j]#distancebetweenatomsiandj129xr-=self.L*round(xr/self.L)#periodicboundaryconditions130r2=xr**2#squaretocomparetocutoff131ifr2self.rc2:#testcutoff132#computeLennard-Jonesinteraction133r2i=1.0/r2134r6i=r2i**3135ff=48.0*r2i*r6i*(r6i-0.5)136#updateforces137self.f[i]+=ff*xr138self.f[j]-=ff*xr139#updateenergy140en+=4.0*r6i*(r6i-1.0)-self.ecut141returnen142143#integrateequationsofmotion144defintegrate(self,t,en):145print(#----Integratingequationsofmotion----)146sumv=0.0147sumv2=0.0148foriinrange(0,self.N):149xx=2.0*self.x[i]-self.xp[i]+self.dt*self.dt*self.f[i]#Verletalgorithm150vi=(xx-self.xp[i])/(2.0*self.dt)#velocity151sumv+=vi#velocitycentreofmass152sumv2+=vi**2#totalkineticenergy153self.xp[i]=self.x[i]#updatepreviouspositions154self.x[i]=xx#updatecurrentpositions155156self.temp=sumv2/(self.dim*self.dN)#instantaneoustemperature157#storeforcalculatingSD158self.sumTemp+=self.temp159self.temps[t]=self.temp160self.etot=(en+0.5*sumv2)/self.dN#totalenergycperparticle161#storeforcalculatingSD162self.sumEtot+=self.etot163self.etots[t]=self.etot164165#printcoordinates166defprintcoords(self,time,cfile):167cfile.write('%d\n'%self.N)168cfile.write('time%10.10f\n'%time)169foriinrange(0,self.N):170cfile.write('C%-8.8f0.00.0\n'%self.x[i])171172173#calculateaverages,etcandprinttofile174defstatistics(self,tfile,efile):175#averages176aveTemp=self.sumTemp/self.dtsteps177aveEtot=self.sumEtot/self.dtsteps178179#standarddeviation180varTemp=0.0181varEtot=0.0182foriinrange(0,self.tsteps):183varTemp+=(self.temps[i]-aveTemp)**2184varEtot+=(self.etots[i]-aveEtot)**2185186sdTemp=math.sqrt(varTemp/self.dtsteps)187sdEtot=math.sqrt(varEtot/self.dtsteps)188tfile.write('#Averagetemperature:%10.2f\n'%aveTemp)189tfile.write('#S