# Copyright 2005-2007 Nanorex, Inc. See LICENSE file for details. """ sim.pyx $Id$ Note: this file is processed by Pyrex to produce sim.c in this directory (not in the build subdirectory). [bruce 060101] Example usage script: make clean; make pyx && python -c "import sim; sim.test()" """ __author__ = "Will" ## import threading #bruce 060101 removed this since I think it's not used now try: import Numeric except ImportError: import numpy as Numeric import unittest cdef extern from "simhelp.c": # note: this produces '#include "simhelp.c"' in generated sim.c file, # but distutils fails to realize there's a dependency on simhelp.c, # so Will added a setup.py dependency to fix that. [bruce 060101] # stuff from globals.c int debug_flags int Iteration int ToMinimize int IterPerFrame int NumFrames int DumpAsText int DumpIntermediateText int PrintFrameNums int OutputFormat int KeyRecordInterval int DirectEvaluate int PrintPotentialEnergy int Interrupted double MinimizeThresholdCutoverRMS double MinimizeThresholdCutoverMax double MinimizeThresholdEndRMS double MinimizeThresholdEndMax double VanDerWaalsCutoffRadius double VanDerWaalsCutoffFactor int EnableElectrostatic int NeighborSearching double ThermostatGamma int UseAMBER int TypeFeedback char *IDKey char *BaseFileName char *InputFileName char *OutputFileName char *TraceFileName char *GromacsOutputBaseName char *PathToCpp char *SystemParametersFileName char *AmberBondedParametersFileName char *AmberNonbondedParametersFileName char *AmberChargesFileName double Dt double Dx double Dmass double Temperature # end of globals.c stuff setWriteTraceCallbackFunc(PyObject) setFrameCallbackFunc(PyObject) getFrame_c() pyrexInitBondTable() void dumpPart() void reinit_globals() everythingElse() everythingDone() cdef char *structCompareHelp() void strcpy(char *, char *) #bruce 051230 guess void reinitSimGlobals(PyObject) specialExceptionIs(PyObject) #void dynamicsMovie_start() #void dynamicsMovie_step() #void dynamicsMovie_finish() void initializeBondTable() double getBondEquilibriumDistance(int element1, int element2, char bondOrder) cdef extern from "string.h": int strcmp(char *s1, char *s2) cdef extern from "stdlib.h": void srand(unsigned int) # wware 060111 a special exception for simulator interruptions class SimulatorInterrupted(Exception): pass specialExceptionIs(SimulatorInterrupted) cdef class _Simulator: """Pyrex permits access to doc strings""" # cdef double *data # bruce 060103 comments: BaseSimulator needs an __init__ method which resets all globals # to their desired initial values, in order to make this correct for successive # uses of one of these objects in one session. ###@@@ # Current code does not support more than one of these objects being active # at one time (trying this will crash); but it's ok to use several in succession # except for the issue of the globals not being reset to their initial values. def __getattr__(self, char *key): if key.startswith('_'): # important optimization (when Python asks for __xxx__) [bruce 060102] raise AttributeError, key if strcmp(key, "debug_flags") == 0: return debug_flags elif strcmp(key, "Iteration") == 0: return Iteration elif strcmp(key, "ToMinimize") == 0: return ToMinimize elif strcmp(key, "IterPerFrame") == 0: return IterPerFrame elif strcmp(key, "NumFrames") == 0: return NumFrames elif strcmp(key, "DumpAsText") == 0: return DumpAsText elif strcmp(key, "DumpIntermediateText") == 0: return DumpIntermediateText elif strcmp(key, "PrintFrameNums") == 0: return PrintFrameNums elif strcmp(key, "OutputFormat") == 0: return OutputFormat elif strcmp(key, "KeyRecordInterval") == 0: return KeyRecordInterval elif strcmp(key, "DirectEvaluate") == 0: return DirectEvaluate elif strcmp(key, "PrintPotentialEnergy") == 0: return PrintPotentialEnergy elif strcmp(key, "Interrupted") == 0: return Interrupted elif strcmp(key, "MinimizeThresholdCutoverRMS") == 0: return MinimizeThresholdCutoverRMS elif strcmp(key, "MinimizeThresholdCutoverMax") == 0: return MinimizeThresholdCutoverMax elif strcmp(key, "MinimizeThresholdEndRMS") == 0: return MinimizeThresholdEndRMS elif strcmp(key, "MinimizeThresholdEndMax") == 0: return MinimizeThresholdEndMax elif strcmp(key, "VanDerWaalsCutoffRadius") == 0: return VanDerWaalsCutoffRadius elif strcmp(key, "VanDerWaalsCutoffFactor") == 0: return VanDerWaalsCutoffFactor elif strcmp(key, "EnableElectrostatic") == 0: return EnableElectrostatic elif strcmp(key, "NeighborSearching") == 0: return NeighborSearching elif strcmp(key, "UseAMBER") == 0: return UseAMBER elif strcmp(key, "TypeFeedback") == 0: return TypeFeedback elif strcmp(key, "ThermostatGamma") == 0: return ThermostatGamma elif strcmp(key, "IDKey") == 0: return IDKey elif strcmp(key, "baseFilename") == 0: #bruce 051230 prevent exception when this is NULL (its default value) if BaseFileName == NULL: # not sure if None would be permitted here # probably it would, but this is better anyway return "" return BaseFileName elif strcmp(key, "InputFileName") == 0: if InputFileName == NULL: # should we raise an AttributeError here? return "" return InputFileName elif strcmp(key, "OutputFileName") == 0: if OutputFileName == NULL: # should we raise an AttributeError here? return "" return OutputFileName elif strcmp(key, "TraceFileName") == 0: if TraceFileName == NULL: # should we raise an AttributeError here? return "" return TraceFileName elif strcmp(key, "GromacsOutputBaseName") == 0: if GromacsOutputBaseName == NULL: # should we raise an AttributeError here? return "" return GromacsOutputBaseName elif strcmp(key, "PathToCpp") == 0: if PathToCpp == NULL: # should we raise an AttributeError here? return "" return PathToCpp elif strcmp(key, "SystemParametersFileName") == 0: if SystemParametersFileName == NULL: # should we raise an AttributeError here? return "" return SystemParametersFileName elif strcmp(key, "AmberBondedParametersFileName") == 0: if AmberBondedParametersFileName == NULL: # should we raise an AttributeError here? return "" return AmberBondedParametersFileName elif strcmp(key, "AmberNonbondedParametersFileName") == 0: if AmberNonbondedParametersFileName == NULL: # should we raise an AttributeError here? return "" return AmberNonbondedParametersFileName elif strcmp(key, "AmberChargesFileName") == 0: if AmberChargesFileName == NULL: # should we raise an AttributeError here? return "" return AmberChargesFileName elif strcmp(key, "Dt") == 0: return Dt elif strcmp(key, "Dx") == 0: return Dx elif strcmp(key, "Dmass") == 0: return Dmass elif strcmp(key, "Temperature") == 0: return Temperature else: raise AttributeError, key def __setattr__(self, char *key, value): if strcmp(key, "debug_flags") == 0: global debug_flags debug_flags = value elif strcmp(key, "Iteration") == 0: global Iteration Iteration = value elif strcmp(key, "ToMinimize") == 0: global ToMinimize ToMinimize = value elif strcmp(key, "IterPerFrame") == 0: global IterPerFrame IterPerFrame = value elif strcmp(key, "NumFrames") == 0: global NumFrames NumFrames = value elif strcmp(key, "DumpAsText") == 0: global DumpAsText DumpAsText = value elif strcmp(key, "DumpIntermediateText") == 0: global DumpIntermediateText DumpIntermediateText = value elif strcmp(key, "PrintFrameNums") == 0: global PrintFrameNums PrintFrameNums = value elif strcmp(key, "OutputFormat") == 0: global OutputFormat OutputFormat = value elif strcmp(key, "KeyRecordInterval") == 0: global KeyRecordInterval KeyRecordInterval = value elif strcmp(key, "DirectEvaluate") == 0: global DirectEvaluate DirectEvaluate = value elif strcmp(key, "PrintPotentialEnergy") == 0: global PrintPotentialEnergy PrintPotentialEnergy = value elif strcmp(key, "Interrupted") == 0: global Interrupted Interrupted = value elif strcmp(key, "MinimizeThresholdCutoverRMS") == 0: global MinimizeThresholdCutoverRMS MinimizeThresholdCutoverRMS = value elif strcmp(key, "MinimizeThresholdCutoverMax") == 0: global MinimizeThresholdCutoverMax MinimizeThresholdCutoverMax = value elif strcmp(key, "MinimizeThresholdEndRMS") == 0: global MinimizeThresholdEndRMS MinimizeThresholdEndRMS = value elif strcmp(key, "MinimizeThresholdEndMax") == 0: global MinimizeThresholdEndMax MinimizeThresholdEndMax = value elif strcmp(key, "VanDerWaalsCutoffRadius") == 0: global VanDerWaalsCutoffRadius VanDerWaalsCutoffRadius = value elif strcmp(key, "VanDerWaalsCutoffFactor") == 0: global VanDerWaalsCutoffFactor VanDerWaalsCutoffFactor = value elif strcmp(key, "EnableElectrostatic") == 0: global EnableElectrostatic EnableElectrostatic = value elif strcmp(key, "NeighborSearching") == 0: global NeighborSearching NeighborSearching = value elif strcmp(key, "UseAMBER") == 0: global UseAMBER UseAMBER = value elif strcmp(key, "TypeFeedback") == 0: global TypeFeedback TypeFeedback = value elif strcmp(key, "ThermostatGamma") == 0: global ThermostatGamma ThermostatGamma = value elif strcmp(key, "IDKey") == 0: global IDKey IDKey = value elif strcmp(key, "baseFilename") == 0: global BaseFileName BaseFileName = value elif strcmp(key, "InputFileName") == 0: global InputFileName InputFileName = value elif strcmp(key, "OutputFileName") == 0: global OutputFileName OutputFileName = value elif strcmp(key, "TraceFileName") == 0: global TraceFileName TraceFileName = value elif strcmp(key, "GromacsOutputBaseName") == 0: global GromacsOutputBaseName GromacsOutputBaseName = value elif strcmp(key, "PathToCpp") == 0: global PathToCpp PathToCpp = value elif strcmp(key, "SystemParametersFileName") == 0: global SystemParametersFileName SystemParametersFileName = value elif strcmp(key, "AmberBondedParametersFileName") == 0: global AmberBondedParametersFileName AmberBondedParametersFileName = value elif strcmp(key, "AmberNonbondedParametersFileName") == 0: global AmberNonbondedParametersFileName AmberNonbondedParametersFileName = value elif strcmp(key, "AmberChargesFileName") == 0: global AmberChargesFileName AmberChargesFileName = value elif strcmp(key, "Dt") == 0: global Dt Dt = value elif strcmp(key, "Dx") == 0: global Dx Dx = value elif strcmp(key, "Dmass") == 0: global Dmass Dmass = value elif strcmp(key, "Temperature") == 0: global Temperature Temperature = value else: raise AttributeError, key def go(self, frame_callback=None, trace_callback=None): "run the simulator loop; optional frame_callback should be None or a callable object" #e if there's an exception in the callback, maybe that should abort the sim run # (requiring change to callback-call-helper in simhelp.c, I think) # and reraise that exception or some other one from this method #bruce 060103 setFrameCallbackFunc(frame_callback) setWriteTraceCallbackFunc(trace_callback) srand(0) everythingElse() # I don't want to bother saving/restoring an old frame_callback, # since I think having a permanent one should be deprecated [bruce 060102] # framebacks are cleared in everythingDone everythingDone() return def structCompare(self): r = structCompare() if r: raise Exception, r def reinitGlobals(self): reinit_globals() def getEquilibriumDistanceForBond(self, element1, element2, order): # element1 and element2 are python ints # order is a python string cdef int int_el1, int_el2 cdef char *c_order int_el1 = element1 int_el2 = element2 c_order = order # initializeBondTable should not require a tracefile, so can be called # at any time. It must have been called before either parsing an mmp # file, or calling getBondEquilibriumDistance(). It checks to see if # it's been called already, so it's cheap if that's the case. It # shouldn't affect the saved warning flags in any way. Warnings are # only printed when bond information is retrieved, and # getBondEquilibriumDistance() avoids the warning code. pyrexInitBondTable() return getBondEquilibriumDistance(int_el1, int_el2, c_order[0]) def getFrame(self): frm = getFrame_c() num_atoms = len(frm) / (3 * 8) array = Numeric.fromstring(frm, Numeric.Float64) return Numeric.resize(array, [num_atoms, 3]) _theSimulator = None def theSimulator(): global _theSimulator if (_theSimulator is None): _theSimulator = _Simulator() return _theSimulator ##################################################################################### # Test code below this point # conventional globals for callback -- they don't have to be used _frameNumber = 0 _frameFreq = 1 _callbackCounter = 0 def _myCallback(last_frame): global _frameNumber _frameNumber = _frameNumber + 1 if (_frameNumber % _frameFreq) == 0: #print "frame %d:" % _frameNumber, getFrame() global _callbackCounter _callbackCounter = _callbackCounter + 1 #bruce 060101 made testsetup a function so it only happens when asked for, # not on every import and sim run (for example, runSim.py doesn't want it) def _testsetup(freq=1): "conventional setup for test functions; returns frame_callback argument for .go method" global _frameNumber, _frameFreq, _callbackCounter _callbackCounter = 0 _frameNumber = 0 _frameFreq = max(1, freq) return _myCallback _tracefile = [ ] def _tracecallback(str): global _tracefile _tracefile.append(str) def _badcallback(str): raise RuntimeError, "This is a bad callback" def _badcallback2(): "This callback should really take an argument" pass ##################################################### def isFileAscii(filename): from string import printable inf = open(filename) R = inf.read(1000) inf.close() nonascii = 0 for x in R: if x not in printable: nonascii = nonascii + 1 return nonascii < 20 class Tests(unittest.TestCase): def setUp(self): global _tracefile _tracefile = [ ] def test_getEquilibriumDistance(self): # try C-C single bond; prints 154.88 assert (theSimulator().getEquilibriumDistanceForBond(6, 6, '1') - 154.88) ** 2 < 0.1 def test_framecallback(self): func = _testsetup(2) m = theSimulator() m.reinitGlobals() m.InputFileName = "tests/minimize/test_h2.mmp" m.OutputFileName = "tests/minimize/test_h2.xyz" m.ToMinimize = 1 m.DumpAsText = 1 m.OutputFormat = 0 m.go(frame_callback=func) assert _callbackCounter == 3, "Callback counter is %d, not 3" %(_callbackCounter) def test_frameAndTraceCallback(self): func = _testsetup(10) d = theSimulator() d.reinitGlobals() d.InputFileName = "tests/rigid_organics/test_C6H10.mmp" d.OutputFileName = "tests/rigid_organics/test_C6H10.dpb" d.ToMinimize = 0 d.DumpAsText = 0 d.OutputFormat = 1 d.go(frame_callback=func, trace_callback=_tracecallback) assert _callbackCounter == 10 def test_traceCallbackWithMotor(self): d = theSimulator() d.reinitGlobals() d.InputFileName = "tests/dynamics/test_0001.mmp" d.OutputFileName = "tests/dynamics/test_0001.dpb" d.ToMinimize = 0 d.DumpAsText = 0 d.OutputFormat = 1 d.go(trace_callback=_tracecallback) # Make sure there is motor information being printed assert len(_tracefile[-5].split()) == 3, \ "Motor info not appearing in trace file:" + _tracefile[18] def test_dpbFileShouldBeBinary(self): d = theSimulator() d.reinitGlobals() d.InputFileName = "tests/dynamics/test_0001.mmp" d.OutputFileName = "tests/dynamics/test_0001.dpb" d.ToMinimize = 0 d.DumpAsText = 0 d.OutputFormat = 1 d.go() # Make sure that the DPB file is really binary assert not isFileAscii("tests/dynamics/test_0001.dpb") ## def test_dynamicsStepStuff(self): ## Dynamics("tests/rigid_organics/test_C6H10.mmp") ## dynamicsMovie_start() ## j = 0 ## for i in range(10000): ## dynamicsMovie_step() ## if (i % 500) == 0: ## #print "Here is frame", i ## #print getFrame() ## j = j + 1 ## dynamicsMovie_finish() ## assert j == 20 ######### Tests that should be expected to fail def test_badCallback1(self): try: d = theSimulator() d.reinitGlobals() d.InputFileName = "tests/dynamics/test_0001.mmp" d.OutputFileName = "tests/dynamics/test_0001.dpb" d.ToMinimize = 0 d.DumpAsText = 0 d.OutputFormat = 1 d.go(trace_callback=_badcallback) assert False, "This test should have failed" except RuntimeError, e: assert e.args[0] == "This is a bad callback" def test_badCallback2(self): try: d = theSimulator() d.reinitGlobals() d.InputFileName = "tests/dynamics/test_0001.mmp" d.OutputFileName = "tests/dynamics/test_0001.dpb" d.ToMinimize = 0 d.DumpAsText = 0 d.OutputFormat = 1 d.go(trace_callback=_badcallback2) assert False, "This test should have failed" except TypeError: pass def test_badCallback3(self): try: d = theSimulator() d.reinitGlobals() d.InputFileName = "tests/dynamics/test_0001.mmp" d.OutputFileName = "tests/dynamics/test_0001.dpb" d.ToMinimize = 0 d.DumpAsText = 0 d.OutputFormat = 1 d.go(trace_callback=42) assert False, "This test should have failed" except RuntimeError, e: assert e.args[0] == "callback is not callable" # # make pyx && python -c "import sim; sim.test()" # def test(): suite = unittest.makeSuite(Tests, 'test') runner = unittest.TextTestRunner() runner.run(suite)