From 0b086051adf1c6dc981cdd7fc723b07ea04907cb Mon Sep 17 00:00:00 2001 From: bannanc <11968541+bannanc@users.noreply.github.com> Date: Wed, 14 Aug 2019 14:40:36 -0600 Subject: [PATCH] add np.array per issue29 --- lib/ase_interface.py | 78 ++++++++++++++++++++++---------------------- 1 file changed, 39 insertions(+), 39 deletions(-) diff --git a/lib/ase_interface.py b/lib/ase_interface.py index 4fc66f9..2f7a467 100644 --- a/lib/ase_interface.py +++ b/lib/ase_interface.py @@ -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) @@ -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 @@ -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)) @@ -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 @@ -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 @@ -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() @@ -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 @@ -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: @@ -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, @@ -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: @@ -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(), @@ -1017,11 +1017,11 @@ 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) @@ -1029,33 +1029,33 @@ def calculate(self, atoms=None, properties=['energy'], 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 = [] @@ -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') @@ -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) @@ -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