pdb2pqr 1.7-2 / usr / share / pdb2pqr / pdb2pka / example.py

"""
    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 `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 xrange(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 xrange(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