-
Notifications
You must be signed in to change notification settings - Fork 0
/
fferror.py
executable file
·83 lines (63 loc) · 2.56 KB
/
fferror.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
#!/usr/bin/env python
import sys
import numpy as np
from numba import jit
import writePDB
from propertyfit.utilities import load_qmfiles
from propertyfit.utilities import angstrom2bohr, hartree2kjmol, number2name
@jit(nopython=True, cache=True)
def charge_esp(rinvmat, testcharges):
natoms = rinvmat.shape[0]
ngridpoints = rinvmat.shape[1]
grid = np.zeros(ngridpoints, dtype=np.float64)
for i in range(natoms):
for j in range(ngridpoints):
grid[j] += testcharges[i] * rinvmat[i,j]
return grid
@jit(nopython=True, cache=True)
def dipole_esp(rinvmat, xyzmat, testdipoles):
natoms = rinvmat.shape[0]
ngridpoints = rinvmat.shape[1]
grid = np.zeros(ngridpoints, dtype=np.float64)
for i in range(natoms):
for j in range(ngridpoints):
for k in range(3):
grid[j] -= testdipoles[i,k] * rinvmat[i,j]**3 * xyzmat[i,j,k]
return grid
if __name__ == "__main__":
qmfile = sys.argv[1]
chargefile = sys.argv[2]
dipolefile = sys.argv[3]
print(qmfile, chargefile, dipolefile)
testcharges = np.loadtxt(chargefile)
testdipoles = np.loadtxt(dipolefile)
structure = load_qmfiles(qmfile)[0]
#PDB out stuff
atom_coords = structure.coordinates
numbers = structure.numbers
labels = [number2name[num] for num in numbers]
print("Computing grid...")
structure.compute_grid(rmin=2.0*angstrom2bohr, rmax=2.0*angstrom2bohr, nsurfaces=1, pointdensity=1.0)
print("Computing rinvmat...")
structure.compute_rinvmat()
print("Computing xyzmat...")
structure.compute_xyzmat()
print("Computing QM ESP...")
structure.compute_qm_esp()
V = np.zeros((structure.grid.shape[0],5))
V[:,0:3] = structure.grid
V[:, 3] = structure.esp_grid_qm
rinvmat = np.copy(structure.rinvmat)
xyzmat = np.copy(structure.xyzmat)
#subtract charge contribution
structure.esp_grid_qm -= charge_esp(rinvmat, testcharges)
V[:,4] = structure.esp_grid_qm
writePDB.writePDB(V, labels, atom_coords, "{}_charge.pdb".format(qmfile))
RMSD = np.sqrt(np.average(structure.esp_grid_qm**2))
print("Charges: RMSD of {} is {} kJ/mol".format(qmfile, RMSD*hartree2kjmol))
#subtract dipole contribution
structure.esp_grid_qm -= dipole_esp(rinvmat, xyzmat, testdipoles)
V[:,4] = structure.esp_grid_qm
writePDB.writePDB(V, labels, atom_coords, "{}_charge+alpha.pdb".format(qmfile))
RMSD = np.sqrt(np.average(structure.esp_grid_qm**2))
print("Charges + Dipoles: RMSD of {} is {} kJ/mol".format(qmfile, RMSD*hartree2kjmol))