1""" Python APBS Driver File
2
3    This module allows a user to run APBS through Python. Use this module if
4    you wish to include APBS in a Python-based application.
5
6    The module mimics the main.c driver that is used in the C version of APBS.
7    The functions which are called are located in apbslib.py, which is
8    automatically generated by SWIG to wrap each APBS function.  See the APBS
9    documentation for more information about each function.
10
11    Todd Dolinsky (todd@ccb.wustl.edu)
12    Nathan Baker (nathan.baker@pnl.gov)
13
14	APBS -- Adaptive Poisson-Boltzmann Solver
15
16	  Nathan A. Baker (nathan.baker@pnl.gov)
17	  Pacific Northwest National Laboratory
18
19	  Additional contributing authors listed in the code documentation.
20
21	Copyright (c) 2010-2011 Battelle Memorial Institute. Developed at the Pacific Northwest National Laboratory, operated by Battelle Memorial Institute, Pacific Northwest Division for the U.S. Department Energy.  Portions Copyright (c) 2002-2010, Washington University in St. Louis.  Portions Copyright (c) 2002-2010, Nathan A. Baker.  Portions Copyright (c) 1999-2002, The Regents of the University of California.  Portions Copyright (c) 1995, Michael Holst
22
23	All rights reserved.
24
25	Redistribution and use in source and binary forms, with or without
26	modification, are permitted provided that the following conditions are met:
27
28	* Redistributions of source code must retain the above copyright notice, this
29	list of conditions and the following disclaimer.
30
31	* Redistributions in binary form must reproduce the above copyright notice,
32	this list of conditions and the following disclaimer in the documentation
33	and/or other materials provided with the distribution.
34
35	* Neither the name of Washington University in St. Louis nor the names of its
36	contributors may be used to endorse or promote products derived from this
37	software without specific prior written permission.
38
39	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
40	"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
41	LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
42	A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
43	CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
44	EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
45	PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
46	PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
47	LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
48	NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
49	SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
50"""
51
52from apbslib import *
53import sys, time
54import string
55from sys import stdout, stderr
56
57__author__ = "Todd Dolinsky, Nathan Baker, Yong Huang"
58__date__ = "September 2009"
59__version__ = "1.3"
60
61Python_kb = 1.3806581e-23
62Python_Na = 6.0221367e+23
63NOSH_MAXMOL = 20
64NOSH_MAXCALC = 20
65
66class APBSError(Exception):
67    """ APBSError class
68
69        The APBSError class inherits off the Exception module and returns
70        a string defining the nature of the error.
71    """
72
73    def __init__(self, value):
74        """
75            Initialize with error message
76
77            Parameters
78                value:  Error Message (string)
79        """
80        self.value = value
81
82    def __str__(self):
83        """
84            Return the error message
85        """
86        return `self.value`
87
88def getHeader():
89    """ Get header information about APBS
90        Returns (header)
91            header: Information about APBS
92    """
93
94    header = "\n\n\
95    ----------------------------------------------------------------------\n\
96    Adaptive Poisson-Boltzmann Solver (APBS)\n\
97    Version 1.3\n\
98    \n\
99    APBS -- Adaptive Poisson-Boltzmann Solver\n\
100    \n\
101    Nathan A. Baker (nathan.baker@pnl.gov)\n\
102    Pacific Northwest National Laboratory\n\
103    \n\
104    Additional contributing authors listed in the code documentation.\n\
105    \n\
106    Copyright (c) 2010-2011 Battelle Memorial Institute. Developed at the Pacific Northwest National Laboratory, operated by Battelle Memorial Institute, Pacific Northwest Division for the U.S. Department Energy.  Portions Copyright (c) 2002-2010, Washington University in St. Louis. Portions Copyright (c) 2002-2010, Nathan A. Baker. Portions Copyright (c) 1999-2002, The Regents of the University of California.  Portions Copyright (c) 1995, Michael Holst\n\
107    \n\
108    All rights reserved.\n\
109    \n\
110    Redistribution and use in source and binary forms, with or without\n\
111    modification, are permitted provided that the following conditions are met:\n\
112    \n\
113    * Redistributions of source code must retain the above copyright notice, this\n\
114      list of conditions and the following disclaimer.\n\
115    \n\
116    * Redistributions in binary form must reproduce the above copyright notice,\n\
117      this list of conditions and the following disclaimer in the documentation\n\
118      and/or other materials provided with the distribution.\n\
119    \n\
120    * Neither the name of Washington University in St. Louis nor the names of its\n\
121      contributors may be used to endorse or promote products derived from this\n\
122      software without specific prior written permission.\n\
123    \n\
124    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS\n\
125    \"AS IS\" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT\n\
126    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR\n\
127    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR\n\
128    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,\n\
129    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,\n\
130    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR\n\
131    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF\n\
132    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING\n\
133    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS\n\
134    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.\n\
135    ----------------------------------------------------------------------\n\
136    \n\n"
137
138    return header
139
140def getUsage():
141    """ Get usage information about running APBS via Python
142        Returns (usage)
143            usage: Text about running APBS via Python
144    """
145
146    usage = "\n\n\
147    ----------------------------------------------------------------------\n\
148    This driver program calculates electrostatic potentials, energies,\n\
149    and forces using both multigrid methods.\n\
150    It is invoked as:\n\n\
151      python main.py apbs.in\n\
152    ----------------------------------------------------------------------\n\n"
153
154    return usage
155
156def main():
157    """ Main driver for testing.  Runs APBS on given input file """
158
159    # Initialize the MALOC library
160    startVio()
161
162    # Initialize variables, arrays
163    com = Vcom_ctor(1)
164    rank = Vcom_rank(com)
165    size = Vcom_size(com)
166    mgparm = MGparm()
167    pbeparm = PBEparm()
168    mem = Vmem_ctor("Main")
169    pbe = new_pbelist(NOSH_MAXMOL)
170    pmg = new_pmglist(NOSH_MAXMOL)
171    pmgp = new_pmgplist(NOSH_MAXMOL)
172    realCenter = double_array(3)
173    totEnergy = []
174    nforce = int_array(NOSH_MAXCALC)
175    atomforce = new_atomforcelist(NOSH_MAXCALC)
176
177    # Start the main timer
178    main_timer_start = time.clock()
179
180    # Check invocation
181    stdout.write(getHeader())
182    if len(sys.argv) != 2:
183        stderr.write("main:  Called with %d arguments!\n" % len(sys.argv))
184        stderr.write(getUsage())
185        raise APBSError, "Incorrect Usage!"
186
187    # Parse the input file
188    nosh = NOsh_ctor(rank, size)
189    input_file = sys.argv[1]
190    stdout.write("Parsing input file %s...\n" % input_file)
191    if NOsh_parseInputFile(nosh, input_file) != 1:
192        stderr.write("main:  Error while parsing input file.\n")
193        raise APBSError, "Error while parsing input file!"
194
195    # Load the molecules using loadMolecules routine
196    # loadMolecule passing NULL as second arg instead of Vparam
197    alist = new_valist(NOSH_MAXMOL)
198    if loadMolecules(nosh,None,alist) != 1:
199        stderr.write("main:  Error while loading molecules. \n")
200        raise APBSError, "Error while loading molecules!"
201
202    # Setup the calculations
203
204    if NOsh_setupElecCalc(nosh, alist) != 1:
205        stderr.write("main: Error while setting up calculations. \n")
206        raise APBSError, "Error while setting up calculations!"
207
208    # Load the necessary maps
209
210    dielXMap = new_gridlist(NOSH_MAXMOL)
211    dielYMap = new_gridlist(NOSH_MAXMOL)
212    dielZMap = new_gridlist(NOSH_MAXMOL)
213    if loadDielMaps(nosh, dielXMap, dielYMap, dielZMap) != 1:
214        stderr.write("Error reading dielectric maps!\n")
215        raise APBSError, "Error reading dielectric maps!"
216
217    kappaMap = new_gridlist(NOSH_MAXMOL)
218    if loadKappaMaps(nosh, kappaMap) != 1:
219        stderr.write("Error reading kappa maps!\n")
220        raise APBSError, "Error reading kappa maps!"
221
222    chargeMap = new_gridlist(NOSH_MAXMOL)
223    if loadChargeMaps(nosh, chargeMap) != 1:
224        stderr.write("Error reading charge maps!\n")
225        raise APBSError, "Error reading charge maps!"
226
227    # Do the calculations
228
229    stdout.write("Preparing to run %d PBE calculations. \n" % nosh.ncalc)
230
231    for icalc in xrange(nosh.ncalc): totEnergy.append(0.0)
232
233    for icalc in xrange(nosh.ncalc):
234        stdout.write("---------------------------------------------\n")
235        calc = NOsh_getCalc(nosh, icalc)
236        mgparm = calc.mgparm
237        pbeparm = calc.pbeparm
238        if calc.calctype != 0:
239            stderr.write("main:  Only multigrid calculations supported!\n")
240            raise APBSError, "Only multigrid calculations supported!"
241
242        for k in range(0, nosh.nelec):
243            if NOsh_elec2calc(nosh,k) >= icalc:
244                break
245
246        name = NOsh_elecname(nosh, k)
247        if name == "":
248            stdout.write("CALCULATION #%d:  MULTIGRID\n" % (icalc+1))
249        else:
250            stdout.write("CALCULATION #%d (%s): MULTIGRID\n" % ((icalc+1),name))
251        stdout.write("Setting up problem...\n")
252
253        # Routine initMG
254
255        if initMG(icalc, nosh, mgparm, pbeparm, realCenter, pbe,
256              alist, dielXMap, dielYMap, dielZMap, kappaMap, chargeMap,
257              pmgp, pmg) != 1:
258            stderr.write("Error setting up MG calculation!\n")
259            raise APBSError, "Error setting up MG calculation!"
260
261        # Print problem parameters if desired (comment out if you want
262        # to minimize output to stdout)
263
264        printMGPARM(mgparm, realCenter)
265        printPBEPARM(pbeparm)
266
267        # Solve the problem : Routine solveMG
268
269        thispmg = get_Vpmg(pmg,icalc)
270
271        if solveMG(nosh, thispmg, mgparm.type) != 1:
272            stderr.write("Error solving PDE! \n")
273            raise APBSError, "Error Solving PDE!"
274
275        # Set partition information : Routine setPartMG
276
277        if setPartMG(nosh, mgparm, thispmg) != 1:
278            stderr.write("Error setting partition info!\n")
279            raise APBSError, "Error setting partition info!"
280
281        # Get the energies - the energy for this calculation
282        # (calculation number icalc) will be stored in the totEnergy array
283
284        ret, totEnergy[icalc] = energyMG(nosh, icalc, thispmg, 0, \
285                                         totEnergy[icalc], 0.0, 0.0, 0.0)
286
287        # Calculate forces
288
289        aforce = get_AtomForce(atomforce, icalc)
290        wrap_forceMG(mem, nosh, pbeparm, mgparm, thispmg, aforce, alist, nforce, icalc)
291
292        # Write out data from MG calculations : Routine writedataMG
293        writedataMG(rank, nosh, pbeparm, thispmg)
294
295        # Write out matrix from MG calculations
296        writematMG(rank, nosh, pbeparm, thispmg)
297
298    # Handle print statements - comment out if limiting output to stdout
299
300    if nosh.nprint > 0:
301        stdout.write("---------------------------------------------\n")
302        stdout.write("PRINT STATEMENTS\n")
303    for iprint in xrange(nosh.nprint):
304        if NOsh_printWhat(nosh, iprint) == NPT_ENERGY:
305            printEnergy(com, nosh, totEnergy, iprint)
306        elif NOsh_printWhat(nosh, iprint) == NPT_FORCE:
307            printForce(com, nosh, nforce, atomforce, iprint)
308        elif NOsh_printWhat(nosh, iprint) == NPT_ELECENERGY:
309            printElecEnergy(com, nosh, totEnergy, iprint)
310        elif NOsh_printWhat(nosh, iprint) == NPT_ELECFORCE:
311            printElecForce(com, nosh, nforce, atomForce, iprint)
312        elif NOsh_printWhat(nosh, iprint) == NPT_APOLENERGY:
313            printApolEnergy(nosh, iprint)
314        elif NOsh_printWhat(nosh, iprint) == NPT_APOLFORCE:
315            printApolForce(com, nosh, nforce, atomForce, iprint)
316        else:
317            stdout.write("Undefined PRINT keyword!\n")
318            break
319
320    stdout.write("----------------------------------------\n")
321    stdout.write("CLEANING UP AND SHUTTING DOWN...\n")
322
323    # Clean up APBS structures
324    killForce(mem, nosh, nforce, atomforce)
325    killEnergy()
326    killMG(nosh, pbe, pmgp, pmg)
327    killChargeMaps(nosh, chargeMap)
328    killKappaMaps(nosh, kappaMap)
329    killDielMaps(nosh, dielXMap, dielYMap, dielZMap)
330    killMolecules(nosh, alist)
331
332    delete_Nosh(nosh)
333
334    # Clean up Python structures
335
336    delete_double_array(realCenter)
337    delete_int_array(nforce)
338    delete_atomforcelist(atomforce)
339    delete_valist(alist)
340    delete_gridlist(dielXMap)
341    delete_gridlist(dielYMap)
342    delete_gridlist(dielZMap)
343    delete_gridlist(kappaMap)
344    delete_gridlist(chargeMap)
345    delete_pmglist(pmg)
346    delete_pmgplist(pmgp)
347    delete_pbelist(pbe)
348
349
350    # Clean up MALOC structures
351    delete_Com(com)
352    delete_Mem(mem)
353    stdout.write("\n")
354    stdout.write("Thanks for using APBS!\n\n")
355
356    # Stop the main timer
357    main_timer_stop = time.clock()
358    stdout.write("Total execution time:  %1.6e sec\n" % (main_timer_stop - main_timer_start))
359
360
361if __name__ == "__main__": main()
362