calulcation on molecular crystals?

[Edward Pyzer-Knapp]

Hi all,

I am writing some code which requires the ability to get energies and forces from a crystal structure (I am doing rigid body minimizations of molecular crystals). I was wondering if someone could help show me how this would look in the OpenMM framework (using the python wrapper would be vastly preferable).
From the manual, I can see that OpenMM can read in box dimensions from a PDB, but I would much rather avoid file I/O.

I was thinking something along the lines of (psuedo-code..)

Code: Select all

system=crystal_system(a,b,c,alpha,beta,gamma,coords)
system.set_forcefield(fooFF)
forces, energy = system.calculate()

I am probably missing something obvious, but would be hugely appreciative of your help.

Ed Pyzer-Knapp

[Lee-Ping Wang]

Hi Ed,

I don't think OpenMM supports general triclinic unit cells at the moment; they have to be orthorhombic (alpha=beta=gamma=90). The function you probably want is simulation.context.setPeriodicBoxVectors(), described here: https://simtk.org/api_docs/openmm/api6_ ... bef702cc6d

Also here's some code for converting lengths/angles to lattice vectors.

Code: Select all

# Container for Bravais lattice vector. Three cell lengths, three angles, three vectors, volume, and TINKER trig functions.
Box = namedtuple('Box',['a','b','c','alpha','beta','gamma','A','B','C','V'])
radian = 180. / np.pi

Box(a,b,c,alpha,beta,gamma,np.array(L1).flatten(),np.array(L2).flatten(),np.array(L3).flatten(),v*a*b*c)

def BuildLatticeFromLengthsAngles(a, b, c, alpha, beta, gamma):
    """ This function takes in three lattice lengths and three lattice angles, and tries to return a complete box specification. """
    alph = alpha*np.pi/180
    bet = beta*np.pi/180
    gamm = gamma*np.pi/180
    v = np.sqrt(1 - cos(alph)**2 - cos(bet)**2 - cos(gamm)**2 + 2*cos(alph)*cos(bet)*cos(gamm))
    Mat = np.matrix([[a, b*cos(gamm), c*cos(bet)],
                  [0, b*sin(gamm), c*((cos(alph)-cos(bet)*cos(gamm))/sin(gamm))],
                  [0, 0, c*v/sin(gamm)]])
    L1 = Mat*np.matrix([[1],[0],[0]])
    L2 = Mat*np.matrix([[0],[1],[0]])
    L3 = Mat*np.matrix([[0],[0],[1]])
    return Box(a,b,c,alpha,beta,gamma,np.array(L1).flatten(),np.array(L2).flatten(),np.array(L3).flatten(),v*a*b*c)