Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

add np.array per issue29 #30

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

add np.array per issue29 #30

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
add np.array per issue29
  • Loading branch information
@bannanc
bannanc committed Aug 14, 2019
commit 0b086051adf1c6dc981cdd7fc723b07ea04907cb
78 changes: 39 additions & 39 deletions lib/ase_interface.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ def calculate(self, atoms=None, properties=['energy'],
#stress_ani = np.zeros((1, 3))
stress_ani = np.zeros(6)


if self.Setup or self.nc.request_setup():
#Setup molecule for MD
natoms = len(self.atoms)
Expand Down Expand Up @@ -241,7 +241,7 @@ def __init__(self, cnstfile, saefile, nnfprefix, Nnet, gpu_list=0, sinet=False,
self.netdict = {'cns': cnstfile,
'sae': saefile,
'nnf': nnfprefix,
'Nn' : self.Nn,
'Nn' : self.Nn,
'epw': enablepairwise}

# number of cores
Expand Down Expand Up @@ -354,7 +354,7 @@ def compute_mean_props(self):
# Compute var sqr
v2 = np.sum(np.sum(np.power(dF, 2), axis=0) / (self.Nn * (self.Nn - 1)))

# Store intermediates
# Store intermediates
self.intermediates = {'var_sqr': v2}

return np.mean(self.E, axis=0), np.mean(self.F, axis=0), np.std(self.E, axis=0) / np.sqrt(float(self.Na)), np.mean(np.std(self.F, axis=0))
Expand Down Expand Up @@ -487,7 +487,7 @@ def compute_mean_props(self):
# Compute var sqr
v2 = np.sum(np.sum(np.power(dF, 2), axis=0) / (self.Nn * (self.Nn - 1)))

# Store intermediates
# Store intermediates
self.intermediates = {'var_sqr': v2}

# Return
Expand Down Expand Up @@ -686,7 +686,7 @@ def calculate(self, atoms=None, properties=['energy'],
# TODO works only for 3D periodic. For 1,2D - update np.zeros((3,3)) part
pbc_inv = (np.linalg.inv(self.atoms.get_cell())).astype(np.float32) if atoms.pbc.all() else np.zeros(
(3, 3), dtype=np.float32)
self.nc.set_cell((self.atoms.get_cell()).astype(np.float32), pbc_inv)
self.nc.set_cell(np.array(self.atoms.get_cell()).astype(np.float32), pbc_inv)

self.Setup = False
## Run this if models are initialized
Expand All @@ -698,7 +698,7 @@ def calculate(self, atoms=None, properties=['energy'],
# TODO works only for 3D periodic. For 1,2D - update np.zeros((3,3)) part
pbc_inv = (np.linalg.inv(self.atoms.get_cell())).astype(np.float32) if atoms.pbc.all() else np.zeros(
(3, 3), dtype=np.float32)
self.nc.set_cell((self.atoms.get_cell()).astype(np.float32), pbc_inv)
self.nc.set_cell(np.array(self.atoms.get_cell()).astype(np.float32), pbc_inv)

## Compute the model properties (you can speed up ASE energy prediction by not doing force backprop unless needed.)
start_time = time.time()
Expand All @@ -721,7 +721,7 @@ def calculate(self, atoms=None, properties=['energy'],
self.results['energy'] = energy
self.Emodel = energy

## If forces are requested store forces
## If forces are requested store forces
if 'forces' in properties:
forces = force

Expand Down Expand Up @@ -751,7 +751,7 @@ def calculate(self, atoms=None, properties=['energy'],
## Compute pairwise if requested
if self.pairwise:
self.add_pairwise(properties)

## Convert energies and forces to requested units
self.results['energy'] = self.energy_conversion * self.results['energy']
if 'forces' in properties:
Expand Down Expand Up @@ -929,10 +929,10 @@ def aniensloader(model, gpu=0, multigpu=False):


# Class of ANI + D3 energies

class ANID3(Calculator):
implemented_properties = ['energy', 'forces', 'stress']

default_parameters = {'xc': 'ani',
'bj': True,
'threebody': True,
Expand All @@ -943,51 +943,51 @@ class ANID3(Calculator):
'rs18': None,
's6': None,
'calculator': None}


nolabel = True

def __init__(self, build=True, gpuid=0, reslist=[], **kwargs):
Calculator.__init__(self, **kwargs)

if build:
anipath = os.path.dirname(__file__)
cnstfile = anipath + '/../ANI-c08f-ntwk/rHCNO-4.6A_16-3.1A_a4-8.params'
saefile = anipath + '/../ANI-c08f-ntwk/sae_6-31gd.dat'
nnfdir = anipath + '/../ANI-c08f-ntwk/networks/'
self.nc = pync.molecule(cnstfile, saefile, nnfdir, gpuid)

self.Setup = True
self.reslist = reslist

def setnc(self, nc):
self.nc = nc

def calculate(self, atoms=None, properties=['energy'],
system_changes=all_changes):
Calculator.calculate(self, atoms, properties, system_changes)

xc = self.parameters.xc.lower()
bj = self.parameters.bj
threebody = self.parameters.threebody
rcut = self.parameters.rcut
rcutcn = self.parameters.rcutcn
calculator = self.parameters.calculator

if bj:
damp = dampbj
else:
damp = damp0

rs6 = s18 = rs18 = s6 = None

try:
rs6, s18, rs18, s6 = damp[xc]
except KeyError:
unknown_functional = True
else:
unknown_functional = False

if self.parameters.s6 is not None:
s6 = self.parameters.s6
if self.parameters.s18 is not None:
Expand All @@ -996,11 +996,11 @@ def calculate(self, atoms=None, properties=['energy'],
rs6 = self.parameters.rs6
if self.parameters.rs18 is not None:
rs18 = self.parameters.rs18

if unknown_functional and None in (s6, s18, rs6, rs18):
raise ValueError("Unknown functional {}! \
Please specify damping parameters.".format(xc))

# D3 calculation part
energy_d3, forces_d3, stress_d3 = d3_calc(
self.atoms.get_atomic_numbers(),
Expand All @@ -1017,45 +1017,45 @@ def calculate(self, atoms=None, properties=['energy'],
pbc=self.atoms.get_pbc(),
bj=bj,
threebody=threebody)

## make up for stress
## TODO
stress_ani = np.zeros(6)

if self.Setup or self.nc.request_setup():
# Setup molecule for MD
natoms = len(self.atoms)
atom_symbols = self.atoms.get_chemical_symbols()
xyz = self.atoms.get_positions()
self.nc.setMolecule(coords=xyz.astype(np.float32), types=atom_symbols)
self.nc.setPBC(self.atoms.get_pbc()[0], self.atoms.get_pbc()[1], self.atoms.get_pbc()[2])

self.Setup = False
else:
xyz = self.atoms.get_positions()
# Set the conformers in NeuroChem
self.nc.setCoordinates(coords=xyz.astype(np.float32))

# TODO works only for 3D periodic. For 1,2D - update np.zeros((3,3)) part
pbc_inv = (np.linalg.inv(self.atoms.get_cell())).astype(np.float32) if atoms.pbc.all() else np.zeros((3, 3),
dtype=np.float32)
self.nc.setCell((self.atoms.get_cell()).astype(np.float32), pbc_inv)
# self.nc.setCell((self.atoms.get_cell()).astype(np.float32),(np.linalg.inv(self.atoms.get_cell())).astype(np.float32))

# start_time2 = time.time()
self.results['energy'] = conv_au_ev * self.nc.energy()[0] + energy_d3
if 'forces' in properties:
forces = -conv_au_ev * self.nc.force() + forces_d3.T

# restrain atoms
for i in self.reslist:
forces[i] = 0.0

self.results['forces'] = forces
self.results['stress'] = conv_au_ev * stress_ani + stress_d3.flat[[0, 4, 8, 5, 2, 1]]
# end_time2 = time.time()
# print('ANI Time:', end_time2 - start_time2)

def __update_neighbors(self):
# print('------------------------')
# szs = []
Expand All @@ -1066,7 +1066,7 @@ def __update_neighbors(self):
# print(an[a])
# szs.append(len(indices))
self.nc.setNeighbors(ind=a, indices=indices.astype(np.int32), offsets=offsets.astype(np.float32))

# indices, offsets = self.nlR.get_neighbors(302)
# f = open('test2.xyz','w')
# f.write(str(len(indices))+'\n')
Expand All @@ -1075,22 +1075,22 @@ def __update_neighbors(self):
# for i, offset in zip(indices, offsets):
# xyz = self.atoms.positions[i]
# f.write(str(an[i]) + ' ' + str(xyz[0]) + ' ' + str(xyz[1]) + ' ' + str(xyz[2]) + '\n')

# print(szs)
# plt.hist(szs, max(szs)-min(szs), normed=1, facecolor='green', alpha=0.75)
# plt.xlabel('Number of neighbors')
# plt.ylabel('Count')
# plt.show()
# print('------------------------')

def get_atomicenergies(self, atoms=None, properties=['energy'],
system_changes=all_changes):
Calculator.calculate(self, atoms, properties, system_changes)

## make up for stress
## TODO
stress_ani = np.zeros((1, 3))

if self.Setup or self.nc.request_setup():
# Setup molecule for MD
natoms = len(self.atoms)
Expand All @@ -1102,9 +1102,9 @@ def get_atomicenergies(self, atoms=None, properties=['energy'],
xyz = self.atoms.get_positions()
# Set the conformers in NeuroChem
self.nc.setCoordinates(coords=xyz.astype(np.float32))

self.nc.energy()

return self.nc.aenergies(True) * conv_au_ev


Expand Down