-
Notifications
You must be signed in to change notification settings - Fork 1
/
MolDigest.py
101 lines (77 loc) · 2.48 KB
/
MolDigest.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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
#Calculate an embeeding for a molecular, such as coulomb matrix
from Mol import *
from Util import *
class MolDigester:
def __init__(self, name_="Coulomb", OType_="Energy"):
self.name = name_
self.OType = OType_
self.lshape = None # output is just the energy
self.eshape = None
def EmbF(self, mol_):
if (self.name =="Coulomb"):
return self.make_cm
else:
raise Exception("Unknown Embedding Function")
return
# def Emb(self, mol_, MakeOutputs=True):
# Ins =
def make_cm(self, mol_):
natoms = mol_.NAtoms()
CM=np.zeros((natoms, natoms))
deri_CM = np.zeros((natoms, natoms, 6))
xyz = (mol_.coords).copy()
ele = (mol_.atoms).copy()
code = """
double dist = 0.0;
double deri[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
for (int j=0; j < natoms; j++) {
for (int k=0; k<natoms; k++) {
if (k==j) {
dist=1.0;
deri[0] =0.0;
deri[1] =0.0;
deri[2] =0.0;
deri[3] =0.0;
deri[4] =0.0;
deri[5] =0.0;
}
else {
dist=sqrt(pow(xyz[j*3+0]-xyz[k*3+0],2) + pow(xyz[j*3+1]-xyz[k*3+1],2) + pow(xyz[j*3+2]-xyz[k*3+2],2));
deri[0] = -(xyz[j*3+0]-xyz[k*3+0])/(dist*dist*dist);
deri[1] = -(xyz[j*3+1]-xyz[k*3+1])/(dist*dist*dist);
deri[2] = -(xyz[j*3+2]-xyz[k*3+2])/(dist*dist*dist);
deri[3] = -deri[0];
deri[4] = -deri[1];
deri[5] = -deri[2];
}
CM[natoms*j+k]= ele[j]*ele[k]/dist;
deri_CM[natoms*j*6+k*6+0]=ele[j]*ele[k]*deri[0];
deri_CM[natoms*j*6+k*6+1]=ele[j]*ele[k]*deri[1];
deri_CM[natoms*j*6+k*6+2]=ele[j]*ele[k]*deri[2];
deri_CM[natoms*j*6+k*6+3]=ele[j]*ele[k]*deri[3];
deri_CM[natoms*j*6+k*6+4]=ele[j]*ele[k]*deri[4];
deri_CM[natoms*j*6+k*6+5]=ele[j]*ele[k]*deri[5];
}
}
"""
res = inline(code, ['CM','natoms','xyz','ele', 'deri_CM'],headers=['<math.h>','<iostream>'], compiler='gcc')
return CM, deri_CM
def GetUpTri(self, CM):
return CM[np.triu_indices(CM.shape[0], 0)]
def EvalDigest(self, mol_):
CM, deri_CM = (self.EmbF(mol_))(mol_)
UpTri = self.GetUpTri(CM)
if self.lshape ==None or self.eshape==None:
self.lshape=1
self.eshape=UpTri.shape[0]
return UpTri, deri_CM
def TrainDigest(self, mol_):
CM, deri_CM = (self.EmbF(mol_))(mol_)
UpTri = self.GetUpTri(CM)
out = mol_.frag_mbe_energy
if self.lshape ==None or self.eshape==None:
self.lshape=1
self.eshape=UpTri.shape[0]
return UpTri, out
def Print(self):
print "Digest name: ", self.name