""" APBS interface for PDB2PQR Todd Dolinsky (todd@ccb.wustl.edu) Washington University in St. Louis Jens Erik Nielsen """ __date__ = "16 August 2005" __author__ = "Todd Dolinsky, Jens Erik Nielsen" import sys import time import string from src import psize from src import inputgen from apbslib import * Python_kb = 1.3806581e-23 Python_Na = 6.0221367e+23 NOSH_MAXMOL = 20 NOSH_MAXCALC = 20 class APBSError(Exception): """ APBSError class The APBSError class inherits off the Exception module and returns a string defining the nature of the error. """ def __init__(self, value): """ Initialize with error message Parameters value: Error Message (string) """ self.value = value def __str__(self): """ Return the error message """ return repr(self.value) def getUnitConversion(): """ Get the unit conversion from kT to kJ/mol Returns factor: The conversion factor (float) """ temp = 298.15 factor = Python_kb/1000.0 * temp * Python_Na return factor def runAPBS(PQR, INPUT): """ Run APBS using PQR and INPUT strings! """ # Initialize the MALOC library startVio() # Initialize variables, arrays com = Vcom_ctor(1) rank = Vcom_rank(com) size = Vcom_size(com) mgparm = MGparm() pbeparm = PBEparm() mem = Vmem_ctor("Main") pbe = new_pbelist(NOSH_MAXMOL) pmg = new_pmglist(NOSH_MAXMOL) pmgp = new_pmgplist(NOSH_MAXMOL) realCenter = double_array(3) totEnergy = [] x = [] y = [] z = [] chg = [] rad = [] # Start the main timer main_timer_start = time.clock() # Parse the input file nosh = NOsh_ctor(rank, size) if not parseInputFromString(nosh, INPUT): stderr.write("main: Error while parsing input file.\n") raise APBSError("Error occurred!") # Load the molecules using Valist_load routine alist = new_valist(NOSH_MAXMOL) #atoms = protein.getAtoms() #protsize = len(atoms) atoms = string.split(PQR,"\n") for i in range(len(atoms)): atom = atoms[i] params = string.split(atom) x.append(float(params[5])) y.append(float(params[6])) z.append(float(params[7])) chg.append(float(params[8])) rad.append(float(params[9])) myAlist = make_Valist(alist,0) Valist_load(myAlist, len(atoms), x,y,z,chg,rad) # Initialize the energy holders for i in range(nosh.ncalc): totEnergy.append(0.0) potList = [] # Initialize the force holders forceList = [] # Load the dieletric maps dielXMap = new_gridlist(NOSH_MAXMOL) dielYMap = new_gridlist(NOSH_MAXMOL) dielZMap = new_gridlist(NOSH_MAXMOL) # Load the kappa maps kappaMap = new_gridlist(NOSH_MAXMOL) # Load the charge maps chargeMap = new_gridlist(NOSH_MAXMOL) # Do the calculations sys.stdout.write("Preparing to run %d PBE calculations. \n" % nosh.ncalc) for icalc in range(nosh.ncalc): sys.stdout.write("---------------------------------------------\n") calc = NOsh_getCalc(nosh, icalc) mgparm = calc.mgparm pbeparm = calc.pbeparm if calc.calctype != 0: sys.stderr.write("main: Only multigrid calculations supported!\n") raise APBSError("Only multigrid calculations supported!") for k in range(0, nosh.nelec): if NOsh_elec2calc(nosh,k) >= icalc: break name = NOsh_elecname(nosh, k+1) if name == "": sys.stdout.write("CALCULATION #%d: MULTIGRID\n" % (icalc+1)) else: sys.stdout.write("CALCULATION #%d (%s): MULTIGRID\n" % ((icalc+1),name)) sys.stdout.write("Setting up problem...\n") # Routine initMG if initMG(icalc, nosh, mgparm, pbeparm, realCenter, pbe, alist, dielXMap, dielYMap, dielZMap, kappaMap, chargeMap, pmgp, pmg) != 1: sys.stderr.write("Error setting up MG calculation!\n") raise APBSError("Error setting up MG calculation!") # Print problem parameters printMGPARM(mgparm, realCenter) printPBEPARM(pbeparm) # Solve the problem : Routine solveMG thispmg = get_Vpmg(pmg,icalc) if solveMG(nosh, thispmg, mgparm.type) != 1: stderr.write("Error solving PDE! \n") raise APBSError("Error Solving PDE!") # Set partition information : Routine setPartMG if setPartMG(nosh, mgparm, thispmg) != 1: sys.stderr.write("Error setting partition info!\n") raise APBSError("Error setting partition info!") ret, totEnergy[icalc] = energyMG(nosh, icalc, thispmg, 0, totEnergy[icalc], 0.0, 0.0, 0.0) # Set partition information # Write out data from MG calculations : Routine writedataMG writedataMG(rank, nosh, pbeparm, thispmg) # Write out matrix from MG calculations writematMG(rank, nosh, pbeparm, thispmg) # GET THE POTENTIALS potentials = getPotentials(nosh, pbeparm, thispmg, myAlist) potList.append(potentials) # Handle print statements if nosh.nprint > 0: sys.stdout.write("---------------------------------------------\n") sys.stdout.write("PRINT STATEMENTS\n") for iprint in range(nosh.nprint): if NOsh_printWhat(nosh, iprint) == NPT_ENERGY: printEnergy(com, nosh, totEnergy, iprint) elif NOsh_printWhat(nosh, iprint) == NPT_FORCE: printForce(com, nosh, nforce, atomforce, iprint) else: sys.stdout.write("Undefined PRINT keyword!\n") break sys.stdout.write("----------------------------------------\n") sys.stdout.write("CLEANING UP AND SHUTTING DOWN...\n") # Clean up APBS structures #killForce(mem, nosh, nforce, atomforce) killEnergy() killMG(nosh, pbe, pmgp, pmg) killChargeMaps(nosh, chargeMap) killKappaMaps(nosh, kappaMap) killDielMaps(nosh, dielXMap, dielYMap, dielZMap) killMolecules(nosh, alist) delete_Nosh(nosh) # Clean up Python structures #ptrfree(nfor) delete_double_array(realCenter) #delete_int_array(nforce) #delete_atomforcelist(atomforce) delete_valist(alist) delete_gridlist(dielXMap) delete_gridlist(dielYMap) delete_gridlist(dielZMap) delete_gridlist(kappaMap) delete_gridlist(chargeMap) delete_pmglist(pmg) delete_pmgplist(pmgp) delete_pbelist(pbe) # Clean up MALOC structures delete_Com(com) delete_Mem(mem) sys.stdout.write("\n") sys.stdout.write("Thanks for using APBS!\n\n") # Stop the main timer main_timer_stop = time.clock() sys.stdout.write("Total execution time: %1.6e sec\n" % (main_timer_stop - main_timer_start)) #Return potentials return potList if __name__ == "__main__": # This would really be something like # atoms = protein.getAtoms() # PQR = protein.printAtoms(atoms) PQR = "ATOM 1 I ION 1 0.000 0.000 0.000 1.00 3.00" size = psize.Psize() size.parseString(PQR) size.setAll() # The actual name doesn't matter since we're loading from the string! input = inputgen.Input("dummy.pqr", size, "mg-auto", 0) # Print out the number of elec statements print("Number of elecs: ", len(input.elecs)) # Let's set the dielectric in the second elec statement input.elecs[1].sdie = 55.55 # Now run APBS potentials = runAPBS(PQR,str(input)) # And print the results! print("Now we have: ", potentials)