# A DRAFT OF A SIMULATION OF 4-ATOM BUTANE # USING THE NEW STRUCTURE from Physical import * from Forces import * from Propagator import * from IO import * from ForceField import * import FTSM io = IO() PHI_DIHEDRAL = 11 PSI_DIHEDRAL = 18 # DEFINE THE NUMBER OF POINTS ON THE STRING numpoints = 8 # PHYSICAL SYSTEMS physarray = [] forcearray = [] proparray = [] temperatures = [] for i in range(0, numpoints): physarray.append([Physical(), Physical()]) forcearray.append([Forces(), Forces()]) temperatures.append([0,0]) for j in range(0, 2): io.readPDBPos(physarray[i][j], "data/alanDipeptideSol/alan_wb10_mineq.pdb") io.readPSF(physarray[i][j], "data/alanDipeptideSol/alan_wb10.psf") io.readPAR(physarray[i][j], "data/alanDipeptideSol/par_all27_prot_lipid.inp") physarray[i][j].bc = "Periodic" physarray[i][j].temperature = 300 physarray[i][j].exclude = "scaled1-4" physarray[i][j].seed = 1234 temperatures[i][j] = physarray[i][j].getTemperature() proparray.append([Propagator(physarray[i][0],forcearray[i][0],io), Propagator(physarray[i][1],forcearray[i][1],io)]) #phys = Physical() #io.readPDBPos(phys, "examples/alanSolStates/alanC7axial_wb5_min_eq.pdb") #io.readPSF(phys, "examples/alanSolStates/alan_wb5.psf") #io.readPAR(phys, "examples/diAlanine/par_all27_prot_lipid.inp") #phys.bc = "Vacuum" #phys.temperature = 300 #phys.exclude = "scaled1-4" #phys.seed = 1234 #forces = Forces() kappa = float(sys.argv[1]) gamma = float(sys.argv[2]) #kappa = 40.0 #gamma = 2000.0 ff = forcearray[0][0].makeForceField(physarray[0][0]) #ff.constrainDihedral(PHI_DIHEDRAL, PSI_DIHEDRAL) ff.bondedForces("badihh") #print ff.forcearray #ff.nonbondedForces("le") ff.nonbondedForces("lc") #ff.setAlgorithm("LennardJonesCoulomb", "SimpleFull") #ff.setSwitching("LennardJones", "C2") #ff.setAlgorithm("LennardJones", "Cutoff") #ff.setCutoff("LennardJones", 20); #ff.setParameters("LennardJones", "switchon=15.0") #ff.setSwitching("CoulombDiElec", "C2") #ff.setAlgorithm("CoulombDiElec", "Cutoff") #ff.setCutoff("CoulombDiElec", 20); #ff.setParameters("CoulombDiElec", "switchon=15.0") #prop = Propagator() # SECOND PHYSICAL SYSTEM #phys2 = Physical() #io.readPDBPos(phys2, "examples/alanSolStates/alanC7axial_wb5_min_eq.pdb") #io.readPSF(phys2, "examples/alanSolStates/alan_wb5.psf") #io.readPAR(phys2, "examples/diAlanine/par_all27_prot_lipid.inp") #phys2.bc = "Vacuum" #phys2.temperature = 300 #phys2.exclude = "scaled1-4" #phys2.seed = 1234 #forces2 = Forces() #prop2 = Propagator() ################################################################### # STEP 1: OBTAIN THE INITIAL STRING S = [] # INITIALIZE TO EMPTY # PROPAGATE ONCE SO OUR CONSTRAINING FORCE CAN GIVE US THE INITIAL # PHI, PSI #phiI = physarray[0][0].phi(PHI_DIHEDRAL) #psiI = physarray[0][0].phi(PSI_DIHEDRAL) phiI = -40*numpy.pi/180 psiI = 130*numpy.pi/180 print phiI, " ", psiI # READ THE SECOND PDB # PROPAGATE ONCE TO GET THE FINAL PHI,PSI #io.readPDBPos(physarray[0][0], "examples/alanSolStates/alanC7equat_wb5_min_eq.pdb") #phiF = physarray[0][0].phi(PHI_DIHEDRAL) #psiF = physarray[0][0].phi(PSI_DIHEDRAL) phiF = -40*numpy.pi/180 psiF = -45*numpy.pi/180 print phiF, " ", psiF # SEPARATE INTO EQUIDISTANT PHI AND PSI dphi = (phiF - phiI) / (numpoints-1) dpsi = (psiF - psiI) / (numpoints-1) S.append([phiI, psiI]) newphi = phiI newpsi = psiI for ii in range(1, numpoints-1): newphi += dphi newpsi += dpsi S.append([newphi, newpsi]) S.append([phiF, psiF]) ff.params['HarmonicDihedral'] = {'kbias':[kappa, kappa], 'dihedralnum':[PHI_DIHEDRAL-1, PSI_DIHEDRAL-1], 'angle':[S[0][0], S[0][1]]} #io.initializePlot('string') #io.pause=1 stringgraph=io.newGraph('Phi', 'Psi') # PRINT INITIAL STRING I0 print "\nI0: ", print S io.plotVector(proparray[0][0],stringgraph,S, rangex=[-numpy.pi, numpy.pi], rangey=[-numpy.pi, numpy.pi]) # SET THE STATE BACK TO THE ORIGINAL PDB #io.readPDBPos(physarray[0][0], "examples/alanSolStates/alanC7axial_wb5_min_eq.pdb") dt = 1.0 for iter in range(0, 200000): # NUMBER OF FTSM ITERATIONS for workpt in range(0, numpoints): # LOOPING OVER POINTS if (iter >= 10000 and iter <= 100000): kappa += (100.-float(sys.argv[1]))/90000. # UPDATE FREE SPACE # USE FIRST SYSTEM TO GET M # USE SECOND SYSTEM TO OBTAIN PHI AND PSI DIFFERENCES # FROM TARGETS M_phiphi = FTSM.M_phi(physarray[workpt][0], PHI_DIHEDRAL) M_psipsi = FTSM.M_phi(physarray[workpt][0], PSI_DIHEDRAL) M_phipsi = FTSM.M_phipsi(physarray[workpt][0], PHI_DIHEDRAL, PSI_DIHEDRAL) #print "DIFFERENCE: ", S[workpt][0]-physarray[workpt][1].phi(PHI_DIHEDRAL), " " , S[workpt][1]-physarray[workpt][1].phi(PSI_DIHEDRAL) S[workpt][0] -= (kappa/gamma)*dt*(M_phiphi*(S[workpt][0]-physarray[workpt][1].phi(PHI_DIHEDRAL)) + M_phipsi*(S[workpt][1] - physarray[workpt][1].phi(PSI_DIHEDRAL))) S[workpt][1] -= (kappa/gamma)*dt*(M_phipsi*(S[workpt][0]-physarray[workpt][1].phi(PHI_DIHEDRAL)) + M_psipsi*(S[workpt][1] - physarray[workpt][1].phi(PSI_DIHEDRAL))) # UPDATE CARTESIAN # Dr. Izaguirre: I have checked and this constraint # is correct. The energy is harmonic, but the force (the gradient) # is not harmonic. In fact it is exactly what is in the paper. proparray[workpt][0].propagate(scheme="NoseHoover", steps=1, dt=dt, forcefield=ff, params={'temp':temperatures[workpt][0], 'inertia':0.01, 'bathpos':0.0}) proparray[workpt][1].propagate(scheme="NoseHoover", steps=1, dt=dt, forcefield=ff, params={'temp':temperatures[workpt][1], 'inertia':0.01, 'bathpos':0.0}) # My own function which sets phi and psi for individual force objects # Saves performance since I only change 'angle', I don't want to define # all new force objects by changing params. if (workpt == numpoints-1): FTSM.setPhiPsi(ff, S[0][0], S[0][1]) else: FTSM.setPhiPsi(ff, S[workpt+1][0], S[workpt+1][1]) #print "\nI"+str(iter+1)+" BEFORE SMOOTH: ",S # REPARAMETERIZE S #S.reverse() #FTSM.switchPhiPsi(S) S = FTSM.arclenChris(FTSM.extractPhi(S), FTSM.extractPsi(S)) #FTSM.switchPhiPsi(S) #S.reverse() print "\nI"+str(iter+1)+": ",S io.plotVector(proparray[0][0],stringgraph,S, rangex=[-numpy.pi, numpy.pi], rangey=[-numpy.pi, numpy.pi])