In [2]:
%load_ext Cython

In [5]:
%%cython --annotate

# Cython force field implementation for trilinear grid

#
# Get all the required declarations
#
include "MMTK/python.pxi"
include "MMTK/numeric.pxi"
include "MMTK/core.pxi"
include "MMTK/universe.pxi"
include 'MMTK/forcefield.pxi'

import gzip

import numpy as np
cimport numpy as np

ctypedef np.float_t float_t
ctypedef np.int_t int_t

#
# The force field term implementation.
# The rules:
#
# - The class must inherit from EnergyTerm.
#
# - EnergyTerm.__init__() must be called with the arguments
#   shown here. The third argument is the name of the EnergyTerm
#   object, the fourth a tuple of the names of all the terms it
#   implements (one object can implement several terms).
#   The assignment to self.eval_func is essential, without it
#   any energy evaluation will crash.
#
# - The function "evaluate" must have exactly the parameter
#   list given in this example.
#
cdef class TrilinearGridTerm(EnergyTerm):
    cdef char* grid_name
    cdef np.ndarray scaling_factor, vals, counts, spacing, hCorner
    cdef int npts, nyz, natoms
    cdef float_t strength, k

    # The __init__ method remembers parameters and loads the potential
    # file. Note that EnergyTerm.__init__ takes care of storing the
    # name and the universe object.
    def __init__(self, universe, spacing, counts, vals, strength,
                 scaling_factor, grid_name):

        EnergyTerm.__init__(self, universe,
                            grid_name, (grid_name,))
        self.eval_func = <void *>TrilinearGridTerm.evaluate

        self.strength = strength
        self.scaling_factor = np.array(scaling_factor, dtype=float)
        self.natoms = len(self.scaling_factor)
        self.grid_name = grid_name

        self.spacing = spacing
        self.counts = counts
        self.vals = vals
        self.nyz = self.counts[1]*self.counts[2]
        self.hCorner = np.array((self.spacing[0]*(self.counts[0]-1),
                        self.spacing[1]*(self.counts[1]-1),
                        self.spacing[2]*(self.counts[2]-1)), dtype=float)
        # To keep atoms within the grid
        self.k = 10000. # kJ/mol nm**2
          
    # This method is called for every single energy evaluation, so make
    # it as efficient as possible. The parameters do_gradients and
    # do_force_constants are flags that indicate if gradients and/or
    # force constants are requested.
    cdef void evaluate(self, PyFFEvaluatorObject *eval,
                       energy_spec *input, energy_data *energy):

        # Input
        cdef vector3 *coordinates
        cdef float_t *scaling_factor
        cdef float_t *vals
        cdef int_t *counts
        cdef float_t *spacing
        cdef float_t *hCorner
        # Output
        cdef float_t gridEnergy
        cdef vector3 *gradients
        # Processing
        cdef int i, ix, iy, iz, atom_index
        cdef float_t vmmm, vmmp, vmpm, vmpp, vpmm, vpmp, vppm, vppp
        cdef float_t vmm, vmp, vpm, vpp, vm, vp
        cdef float_t fx, fy, fz, ax, ay, az
        cdef float_t dvdx, dvdy, dvdz

        gridEnergy = 0
        coordinates = <vector3 *>input.coordinates.data

        # Pointers to numpy arrays for faster indexing
        scaling_factor = <float_t *>self.scaling_factor.data
        vals = <float_t *>self.vals.data
        counts = <int_t *>self.counts.data
        spacing = <float_t *>self.spacing.data
        hCorner = <float_t *>self.hCorner.data

        # Initialize variables
        if energy.gradients != NULL:
          gradients = <vector3 *>(<PyArrayObject *> energy.gradients).data
      
        indicies = [ai for ai in range(self.natoms) if self.scaling_factor[ai]!=0]
        for atom_index in indicies:
          # Check to make sure coordinate is in grid
          if (coordinates[atom_index][0]>0 and 
              coordinates[atom_index][1]>0 and 
              coordinates[atom_index][2]>0 and
              coordinates[atom_index][0]<hCorner[0] and
              coordinates[atom_index][1]<hCorner[1] and
              coordinates[atom_index][2]<hCorner[2]):

            # Index within the grid
            ix = int(coordinates[atom_index][0]/spacing[0])
            iy = int(coordinates[atom_index][1]/spacing[1])
            iz = int(coordinates[atom_index][2]/spacing[2])
            
            i = ix*self.nyz + iy*counts[2] + iz

            # Corners of the box surrounding the point
            vmmm = vals[i]
            vmmp = vals[i+1]
            vmpm = vals[i+counts[2]]
            vmpp = vals[i+counts[2]+1]

            vpmm = vals[i+self.nyz]
            vpmp = vals[i+self.nyz+1]
            vppm = vals[i+self.nyz+counts[2]]
            vppp = vals[i+self.nyz+counts[2]+1]
            
            # Fraction within the box
            fx = (coordinates[atom_index][0] - (ix*spacing[0]))/spacing[0]
            fy = (coordinates[atom_index][1] - (iy*spacing[1]))/spacing[1]
            fz = (coordinates[atom_index][2] - (iz*spacing[2]))/spacing[2]
            
            # Fraction ahead
            ax = 1 - fx
            ay = 1 - fy
            az = 1 - fz
      
            # Trilinear interpolation for energy
            vmm = az*vmmm + fz*vmmp
            vmp = az*vmpm + fz*vmpp
            vpm = az*vpmm + fz*vpmp
            vpp = az*vppm + fz*vppp
            
            vm = ay*vmm + fy*vmp
            vp = ay*vpm + fy*vpp
            
            gridEnergy += scaling_factor[atom_index]*(ax*vm + fx*vp)
          
            if energy.gradients != NULL:
              # x coordinate
              dvdx = -vm + vp
              # y coordinate
              dvdy = (-vmm + vmp)*ax + (-vpm + vpp)*fx
              # z coordinate
              dvdz = ((-vmmm + vmmp)*ay + (-vmpm + vmpp)*fy)*ax + ((-vpmm + vpmp)*ay + (-vppm + vppp)*fy)*fx
            
              gradients[atom_index][0] += self.strength*scaling_factor[atom_index]*dvdx/spacing[0]
              gradients[atom_index][1] += self.strength*scaling_factor[atom_index]*dvdy/spacing[1]
              gradients[atom_index][2] += self.strength*scaling_factor[atom_index]*dvdz/spacing[2]
          else:
            for i in range(3):
              if (coordinates[atom_index][i]<0):
                gridEnergy += self.k*coordinates[atom_index][i]**2/2.
                if energy.gradients != NULL:
                  gradients[atom_index][i] += self.k*coordinates[atom_index][i]
              elif (coordinates[atom_index][i]>hCorner[i]):
                gridEnergy += self.k*(coordinates[atom_index][i]-hCorner[i])**2/2.
                if energy.gradients != NULL:
                  gradients[atom_index][i] += self.k*(coordinates[atom_index][i]-hCorner[i])

        energy.energy_terms[self.index] = gridEnergy*self.strength


warning: /Users/dminh/Installers/MMTK-2.7.9/Include/MMTK/numeric.pxi:22:12: Non-trivial type declarators in shared declaration (e.g. mix of pointers and values). Each pointer declaration should be on its own line.
warning: /Users/dminh/Installers/MMTK-2.7.9/Include/MMTK/numeric.pxi:22:25: Non-trivial type declarators in shared declaration (e.g. mix of pointers and values). Each pointer declaration should be on its own line.
warning: /Users/dminh/Installers/MMTK-2.7.9/Include/MMTK/numeric.pxi:30:17: Non-trivial type declarators in shared declaration (e.g. mix of pointers and values). Each pointer declaration should be on its own line.
warning: /Users/dminh/Installers/MMTK-2.7.9/Include/MMTK/numeric.pxi:30:30: Non-trivial type declarators in shared declaration (e.g. mix of pointers and values). Each pointer declaration should be on its own line.
---------------------------------------------------------------------------
CompileError                              Traceback (most recent call last)
<ipython-input-5-09fc42960041> in <module>()
----> 1 get_ipython().run_cell_magic(u'cython', u'--annotate', u'\n# Cython force field implementation for trilinear grid\n\n#\n# Get all the required declarations\n#\ninclude "/Users/dminh/Installers/MMTK-2.7.9/Include/MMTK/python.pxi"\ninclude "/Users/dminh/Installers/MMTK-2.7.9/Include/MMTK/numeric.pxi"\ninclude "/Users/dminh/Installers/MMTK-2.7.9/Include/MMTK/core.pxi"\ninclude "/Users/dminh/Installers/MMTK-2.7.9/Include/MMTK/universe.pxi"\ninclude \'/Users/dminh/Installers/MMTK-2.7.9/Include/MMTK/forcefield.pxi\'\n\nimport gzip\n\nimport numpy as np\ncimport numpy as np\n\nctypedef np.float_t float_t\nctypedef np.int_t int_t\n\n#\n# The force field term implementation.\n# The rules:\n#\n# - The class must inherit from EnergyTerm.\n#\n# - EnergyTerm.__init__() must be called with the arguments\n#   shown here. The third argument is the name of the EnergyTerm\n#   object, the fourth a tuple of the names of all the terms it\n#   implements (one object can implement several terms).\n#   The assignment to self.eval_func is essential, without it\n#   any energy evaluation will crash.\n#\n# - The function "evaluate" must have exactly the parameter\n#   list given in this example.\n#\ncdef class TrilinearGridTerm(EnergyTerm):\n    cdef char* grid_name\n    cdef np.ndarray scaling_factor, vals, counts, spacing, hCorner\n    cdef int npts, nyz, natoms\n    cdef float_t strength, k\n\n    # The __init__ method remembers parameters and loads the potential\n    # file. Note that EnergyTerm.__init__ takes care of storing the\n    # name and the universe object.\n    def __init__(self, universe, spacing, counts, vals, strength,\n                 scaling_factor, grid_name):\n\n        EnergyTerm.__init__(self, universe,\n                            grid_name, (grid_name,))\n        self.eval_func = <void *>TrilinearGridTerm.evaluate\n\n        self.strength = strength\n        self.scaling_factor = np.array(scaling_factor, dtype=float)\n        self.natoms = len(self.scaling_factor)\n        self.grid_name = grid_name\n\n        self.spacing = spacing\n        self.counts = counts\n        self.vals = vals\n        self.nyz = self.counts[1]*self.counts[2]\n        self.hCorner = np.array((self.spacing[0]*(self.counts[0]-1),\n                        self.spacing[1]*(self.counts[1]-1),\n                        self.spacing[2]*(self.counts[2]-1)), dtype=float)\n        # To keep atoms within the grid\n        self.k = 10000. # kJ/mol nm**2\n          \n    # This method is called for every single energy evaluation, so make\n    # it as efficient as possible. The parameters do_gradients and\n    # do_force_constants are flags that indicate if gradients and/or\n    # force constants are requested.\n    cdef void evaluate(self, PyFFEvaluatorObject *eval,\n                       energy_spec *input, energy_data *energy):\n\n        # Input\n        cdef vector3 *coordinates\n        cdef float_t *scaling_factor\n        cdef float_t *vals\n        cdef int_t *counts\n        cdef float_t *spacing\n        cdef float_t *hCorner\n        # Output\n        cdef float_t gridEnergy\n        cdef vector3 *gradients\n        # Processing\n        cdef int i, ix, iy, iz, atom_index\n        cdef float_t vmmm, vmmp, vmpm, vmpp, vpmm, vpmp, vppm, vppp\n        cdef float_t vmm, vmp, vpm, vpp, vm, vp\n        cdef float_t fx, fy, fz, ax, ay, az\n        cdef float_t dvdx, dvdy, dvdz\n\n        gridEnergy = 0\n        coordinates = <vector3 *>input.coordinates.data\n\n        # Pointers to numpy arrays for faster indexing\n        scaling_factor = <float_t *>self.scaling_factor.data\n        vals = <float_t *>self.vals.data\n        counts = <int_t *>self.counts.data\n        spacing = <float_t *>self.spacing.data\n        hCorner = <float_t *>self.hCorner.data\n\n        # Initialize variables\n        if energy.gradients != NULL:\n          gradients = <vector3 *>(<PyArrayObject *> energy.gradients).data\n      \n        indicies = [ai for ai in range(self.natoms) if self.scaling_factor[ai]!=0]\n        for atom_index in indicies:\n          # Check to make sure coordinate is in grid\n          if (coordinates[atom_index][0]>0 and \n              coordinates[atom_index][1]>0 and \n              coordinates[atom_index][2]>0 and\n              coordinates[atom_index][0]<hCorner[0] and\n              coordinates[atom_index][1]<hCorner[1] and\n              coordinates[atom_index][2]<hCorner[2]):\n\n            # Index within the grid\n            ix = int(coordinates[atom_index][0]/spacing[0])\n            iy = int(coordinates[atom_index][1]/spacing[1])\n            iz = int(coordinates[atom_index][2]/spacing[2])\n            \n            i = ix*self.nyz + iy*counts[2] + iz\n\n            # Corners of the box surrounding the point\n            vmmm = vals[i]\n            vmmp = vals[i+1]\n            vmpm = vals[i+counts[2]]\n            vmpp = vals[i+counts[2]+1]\n\n            vpmm = vals[i+self.nyz]\n            vpmp = vals[i+self.nyz+1]\n            vppm = vals[i+self.nyz+counts[2]]\n            vppp = vals[i+self.nyz+counts[2]+1]\n            \n            # Fraction within the box\n            fx = (coordinates[atom_index][0] - (ix*spacing[0]))/spacing[0]\n            fy = (coordinates[atom_index][1] - (iy*spacing[1]))/spacing[1]\n            fz = (coordinates[atom_index][2] - (iz*spacing[2]))/spacing[2]\n            \n            # Fraction ahead\n            ax = 1 - fx\n            ay = 1 - fy\n            az = 1 - fz\n      \n            # Trilinear interpolation for energy\n            vmm = az*vmmm + fz*vmmp\n            vmp = az*vmpm + fz*vmpp\n            vpm = az*vpmm + fz*vpmp\n            vpp = az*vppm + fz*vppp\n            \n            vm = ay*vmm + fy*vmp\n            vp = ay*vpm + fy*vpp\n            \n            gridEnergy += scaling_factor[atom_index]*(ax*vm + fx*vp)\n          \n            if energy.gradients != NULL:\n              # x coordinate\n              dvdx = -vm + vp\n              # y coordinate\n              dvdy = (-vmm + vmp)*ax + (-vpm + vpp)*fx\n              # z coordinate\n              dvdz = ((-vmmm + vmmp)*ay + (-vmpm + vmpp)*fy)*ax + ((-vpmm + vpmp)*ay + (-vppm + vppp)*fy)*fx\n            \n              gradients[atom_index][0] += self.strength*scaling_factor[atom_index]*dvdx/spacing[0]\n              gradients[atom_index][1] += self.strength*scaling_factor[atom_index]*dvdy/spacing[1]\n              gradients[atom_index][2] += self.strength*scaling_factor[atom_index]*dvdz/spacing[2]\n          else:\n            for i in range(3):\n              if (coordinates[atom_index][i]<0):\n                gridEnergy += self.k*coordinates[atom_index][i]**2/2.\n                if energy.gradients != NULL:\n                  gradients[atom_index][i] += self.k*coordinates[atom_index][i]\n              elif (coordinates[atom_index][i]>hCorner[i]):\n                gridEnergy += self.k*(coordinates[atom_index][i]-hCorner[i])**2/2.\n                if energy.gradients != NULL:\n                  gradients[atom_index][i] += self.k*(coordinates[atom_index][i]-hCorner[i])\n\n        energy.energy_terms[self.index] = gridEnergy*self.strength')

/Users/dminh/Applications/miniconda2/envs/algdock/lib/python2.7/site-packages/IPython/core/interactiveshell.pyc in run_cell_magic(self, magic_name, line, cell)
   2118             magic_arg_s = self.var_expand(line, stack_depth)
   2119             with self.builtin_trap:
-> 2120                 result = fn(magic_arg_s, cell)
   2121             return result
   2122 

<decorator-gen-126> in cython(self, line, cell)

/Users/dminh/Applications/miniconda2/envs/algdock/lib/python2.7/site-packages/IPython/core/magic.pyc in <lambda>(f, *a, **k)
    191     # but it's overkill for just that one bit of state.
    192     def magic_deco(arg):
--> 193         call = lambda f, *a, **k: f(*a, **k)
    194 
    195         if callable(arg):

/Users/dminh/Applications/miniconda2/envs/algdock/lib/python2.7/site-packages/Cython/Build/IpythonMagic.pyc in cython(self, line, cell)
    276             build_extension.build_temp = os.path.dirname(pyx_file)
    277             build_extension.build_lib  = lib_dir
--> 278             build_extension.run()
    279             self._code_cache[key] = module_name
    280 

/Users/dminh/Applications/miniconda2/envs/algdock/lib/python2.7/distutils/command/build_ext.pyc in run(self)
    337 
    338         # Now actually compile and link everything.
--> 339         self.build_extensions()
    340 
    341     def check_extensions_list(self, extensions):

/Users/dminh/Applications/miniconda2/envs/algdock/lib/python2.7/distutils/command/build_ext.pyc in build_extensions(self)
    446 
    447         for ext in self.extensions:
--> 448             self.build_extension(ext)
    449 
    450     def build_extension(self, ext):

/Users/dminh/Applications/miniconda2/envs/algdock/lib/python2.7/distutils/command/build_ext.pyc in build_extension(self, ext)
    496                                          debug=self.debug,
    497                                          extra_postargs=extra_args,
--> 498                                          depends=ext.depends)
    499 
    500         # XXX -- this is a Vile HACK!

/Users/dminh/Applications/miniconda2/envs/algdock/lib/python2.7/distutils/ccompiler.pyc in compile(self, sources, output_dir, macros, include_dirs, debug, extra_preargs, extra_postargs, depends)
    572             except KeyError:
    573                 continue
--> 574             self._compile(obj, src, ext, cc_args, extra_postargs, pp_opts)
    575 
    576         # Return *all* object filenames, not just the ones we just built.

/Users/dminh/Applications/miniconda2/envs/algdock/lib/python2.7/distutils/unixccompiler.pyc in _compile(self, obj, src, ext, cc_args, extra_postargs, pp_opts)
    120                        extra_postargs)
    121         except DistutilsExecError, msg:
--> 122             raise CompileError, msg
    123 
    124     def create_static_lib(self, objects, output_libname,

CompileError: command 'gcc' failed with exit status 1

In [1]:
import AlGDock

from MMTK import *
import Interpolation
from MMTK.ForceFields.ForceFieldTest import gradientTest

universe = InfiniteUniverse()

universe.atom1 = Atom('C', position=Vector(1.1, 0.5, 1.5))
universe.atom1.test_charge = 1.
universe.atom2 = Atom('C', position=Vector(1.553, 1.724, 1.464))
universe.atom2.test_charge = -0.2

param_sets = [\
  {'interpolation_type':'Trilinear', 'inv_power':None, 'energy_thresh':-1.0}]

steps = 50000
import numpy as np

from collections import OrderedDict
Es = OrderedDict()
x=np.linspace(1.35,1.6,steps)

for params in param_sets:
  print
  print params
  print
  
  ForceField = Interpolation.InterpolationForceField(\
    '../../../Example/grids/LJa.nc',
    interpolation_type=params['interpolation_type'],
    inv_power=params['inv_power'],
    energy_thresh=params['energy_thresh'],
    scaling_property='test_charge')

  universe.setForceField(ForceField)

  universe.atom1.setPosition(Vector(x[0],0.5,1.5))

  print 'Energy Terms:'
  print universe.energyTerms()
  e, g = universe.energyAndGradients()
  print 'Gradient on Atom 1'
  print g[universe.atom1]
  print 'Gradient on Atom 2'
  print g[universe.atom2]

  print 'Gradient Test'
  gradientTest(universe)

  import time
  start_time = time.time()
  
  key = params['interpolation_type']
  if params['inv_power'] is not None:
    key += ', x**%d'%params['inv_power']
  if params['energy_thresh']>0:
    key += ', e<%f'%params['energy_thresh']

  Es[key] = np.zeros((steps,1))
  for n in range(steps):
    universe.atom1.setPosition(Vector(x[n],0.5,1.5))
    e, g = universe.energyAndGradients()
    Es[key][n] = e
  print 'Time to do %d energy and gradient evaluations: %f s'%(\
    steps, time.time()-start_time)

%matplotlib inline
import matplotlib.pyplot as plt
for key in Es.keys():
  plt.plot(x,Es[key])
plt.legend(Es.keys())


{'interpolation_type': 'Trilinear', 'energy_thresh': -1.0, 'inv_power': None}

Energy Terms:
{'Interpolation': -0.5035438268671999}
Gradient on Atom 1
[2.4492360000000035, 4.4652519999999996, 2.356524000000002]
Gradient on Atom 2
[-1.1748288704000007, 0.5334514112000035, -0.6846025087999981]
Gradient Test
Energy:  -0.503543826867
Atom carbon
[2.4492360000000035, 4.4652519999999996, 2.356524000000002]
[3.106016000001266, 6.027314000000561, 2.9514220000015357]
Atom carbon
[-1.1748288704000007, 0.5334514112000035, -0.6846025087999981]
[-1.174828870399991, 0.53345141119987, -0.6846025088003138]
Time to do 50000 energy and gradient evaluations: 1.117123 s
Out[1]:
<matplotlib.legend.Legend at 0x104e111d0>

In [ ]: