From acf2c65f3011e3ffc14510982b4a7ceb17e08b70 Mon Sep 17 00:00:00 2001 From: qzhu2017 Date: Sat, 2 Mar 2024 12:12:46 -0500 Subject: [PATCH] fix the issue of representation.from_string --- pyxtal/interface/dftb.py | 48 +++++++++++++++++------------------ pyxtal/interface/lammpslib.py | 32 +++++++++++------------ pyxtal/representation.py | 30 ++++++++++++++-------- 3 files changed, 59 insertions(+), 51 deletions(-) diff --git a/pyxtal/interface/dftb.py b/pyxtal/interface/dftb.py index ad704789..6582862c 100644 --- a/pyxtal/interface/dftb.py +++ b/pyxtal/interface/dftb.py @@ -73,7 +73,7 @@ def make_Hamiltonian(skf_dir, atom_types, disp, kpts, write_band=False, use_omp= 'Analysis_MullikenAnalysis': 'No', 'Analysis_CalculateForces': 'Yes', } - if write_band: + if write_band: kwargs['Analysis_WriteBandOut'] = 'Yes' if use_omp: kwargs['Parallel_'] = '' @@ -81,9 +81,9 @@ def make_Hamiltonian(skf_dir, atom_types, disp, kpts, write_band=False, use_omp= if dispersion is not None: kwargs['Hamiltonian_Dispersion'] = dispersion - if skf_dir.find('3ob') > 0: + if skf_dir.find('3ob') > 0: calc_type = '3ob' - elif skf_dir.find('mio') > 0: + elif skf_dir.find('mio') > 0: calc_type = 'mio' elif skf_dir.find('pbc') > 0: calc_type = 'pbc' @@ -97,16 +97,16 @@ def make_Hamiltonian(skf_dir, atom_types, disp, kpts, write_band=False, use_omp= HD = {"Br": -0.0573, "C": -0.1492, "N": -0.1535, - "Ca": -0.0340, + "Ca": -0.0340, "Na": -0.0454, - "Cl": -0.0697, + "Cl": -0.0697, "Zn": -0.03, "O": -0.1575, "F": -0.1623, "P": -0.14, - "H": -0.1857, + "H": -0.1857, "S": -0.11, - "I": -0.0433, + "I": -0.0433, "K": -0.0339, } strs = '{' @@ -144,7 +144,7 @@ def make_Hamiltonian(skf_dir, atom_types, disp, kpts, write_band=False, use_omp= kwargs['Hamiltonian_MaxAngularMomentum_'+ele]='d' else: raise RuntimeError(calc_type, "doesnot support", ele) - + #DFTB2 #pbc-0-3 @@ -183,7 +183,7 @@ def DFTB_relax(struc, skf_dir, opt_cell=False, step=500, \ struc.set_calculator(calc) # impose symmetry - if symmetrize: struc.set_constraint(FixSymmetry(struc)) + if symmetrize: struc.set_constraint(FixSymmetry(struc)) # impose cell constraints if opt_cell: @@ -210,8 +210,8 @@ def DFTB_SCF(struc, skf_dir, kresol=0.10, folder='tmp', disp=None, filename=None Args: struc: ase atoms object - skf_dir: - kresol: + skf_dir: + kresol: filename: band structure output """ if not os.path.exists(folder): @@ -257,10 +257,10 @@ class DFTB(): use_omp: True/False """ - def __init__(self, struc, skf_dir, - disp = 'D3', - kresol = 0.10, - folder = 'tmp', + def __init__(self, struc, skf_dir, + disp = 'D3', + kresol = 0.10, + folder = 'tmp', label = 'test', prefix = 'geo_final', use_omp = False, @@ -276,7 +276,7 @@ def __init__(self, struc, skf_dir, self.prefix = prefix self.use_omp = use_omp if not os.path.exists(self.folder): - os.makedirs(self.folder) + os.makedirs(self.folder) def get_calculator(self, mode, step=500, ftol=1e-3, FixAngles=False, eVperA=True, md_params={}): """ @@ -295,7 +295,7 @@ def get_calculator(self, mode, step=500, ftol=1e-3, FixAngles=False, eVperA=True if eVperA: ftol *= 0.194469064593167E-01 atom_types = set(self.struc.get_chemical_symbols()) kwargs = make_Hamiltonian(self.skf_dir, atom_types, self.disp, self.kpts, use_omp=self.use_omp) - + if mode in ['relax', 'vc-relax']: #kwargs['Driver_'] = 'ConjugateGradient' kwargs['Driver_'] = 'GeometryOptimization' @@ -305,7 +305,7 @@ def get_calculator(self, mode, step=500, ftol=1e-3, FixAngles=False, eVperA=True kwargs['Driver_OutputPrefix'] = self.prefix kwargs['Driver_Convergence_'] = '' kwargs['Driver_Convergence_GradElem'] = ftol - + if mode == 'vc-relax': kwargs['Driver_MovedAtoms'] = "1:-1" kwargs['Driver_LatticeOpt'] = "Yes" @@ -343,7 +343,7 @@ def get_calculator(self, mode, step=500, ftol=1e-3, FixAngles=False, eVperA=True kwargs['Driver_Barostat_Pressure [Pa]'] = dicts['pressure'] kwargs['Driver_Barostat_Timescale [ps]'] = 0.1 - + calc = Dftb(label=self.label, #run_manyDftb_steps=True, atoms=self.struc, @@ -383,10 +383,10 @@ def run(self, mode, step=500, ftol=1e-3, FixAngles=False, md_params={}): class Dftb(FileIOCalculator): """ This module defines a FileIOCalculator for DFTB+ - + http://www.dftbplus.org/ http://www.dftb.org/ - + Initial development: markus.kaukonen@iki.fi Modified by QZ to avoid the I/O load """ @@ -654,7 +654,7 @@ def read_forces(self): # reaches outside of the for loop index_force_begin = -1 index_force_end = -1 - + # Force line indexes fstring = 'forces ' for iline, line in enumerate(self.lines): @@ -670,8 +670,8 @@ def read_forces(self): gradients.append([float(word[k]) for k in range(0, 3)]) gradients = np.array(gradients)* Hartree / Bohr - return gradients - + return gradients + def read_eigenvalues(self): """ Read Eigenvalues from dftb output file (results.tag). Unfortunately, the order seems to be scrambled. """ diff --git a/pyxtal/interface/lammpslib.py b/pyxtal/interface/lammpslib.py index 46daea78..40922afa 100755 --- a/pyxtal/interface/lammpslib.py +++ b/pyxtal/interface/lammpslib.py @@ -32,7 +32,7 @@ def opt_lammpslib(struc, lmp, parameters=None, mask=None, logfile='-', #check_symmetry(si, 1.0e-6, verbose=True) struc.set_calculator(lammps) - struc.set_constraint(FixSymmetry(struc)) + struc.set_constraint(FixSymmetry(struc)) if opt_cell: ecf = ExpCellFilter(struc, mask) if method == 'FIRE': @@ -84,17 +84,17 @@ def run_lammpslib(struc, lmp, parameters=None, path='tmp', \ f"min_style {min_style}", f"minimize 1e-15 1e-15 {steps} {steps}", ] - + else: parameter0 += ['run 0'] - + lammps = LAMMPSlib(lmp=lmp, lmpcmds=parameter0, molecule=molecule, \ lmp_file=lmp_file, log_file='lammps.log', path=path) - + struc.set_calculator(lammps) - + Eng = struc.get_potential_energy() - + return struc, Eng @@ -159,7 +159,7 @@ def __init__(self, lmp, lmpcmds=None, path='tmp', molecule=False, if not self.molecule and para.find('pair_coeff') >= 0: tmp = para.split() filename = tmp[3] - filename1 = self.ffpath + '/' + filename + filename1 = self.ffpath + '/' + filename para = para.replace(filename, filename1) self.paras.append(para) @@ -172,7 +172,7 @@ def calculate(self, atoms): boundary = '' for i in range(3): if atoms.pbc[i]: - boundary += 'p ' + boundary += 'p ' else: boundary += 'f ' if boundary in ['f f p ', 'p p f ']: #needs some work later @@ -203,8 +203,8 @@ def calculate(self, atoms): pos = np.array( [x for x in self.lmp.gather_atoms("x", 1, 3)]).reshape(-1, 3) - - self.energy = self.lmp.extract_variable('pe', None, 0) + + self.energy = self.lmp.extract_variable('pe', None, 0) #print('update lammps energy') xlo = self.lmp.extract_global("boxxlo", 1) @@ -274,7 +274,7 @@ def write_lammps_in(self): fh.write('boundary p p s\n') else: fh.write('boundary {:s}\n'.format(self.boundary)) - fh.write('atom_modify sort 0 0.0\n') + fh.write('atom_modify sort 0 0.0\n') if not self.molecule: fh.write('read_data {:s}\n'.format(self.lammps_data)) fh.write('\n### interactions\n') @@ -297,7 +297,7 @@ def write_lammps_in(self): tmp += LAMMPS_collections().dict['GB_opt'] for i in tmp: fh.write(i) - + if self.calc_type.find('GB') >= 0: fh.write('thermo_style custom pe pxx pyy pzz pyz pxz pxy c_eatoms\n') else: @@ -306,7 +306,7 @@ def write_lammps_in(self): fh.write('thermo 1\n') fh.write('run 1\n') - #fh.write('print "__end_of_ase_invoked_calculation__"\n') + #fh.write('print "__end_of_ase_invoked_calculation__"\n') def write_lammps_data(self, atoms): atom_types = np.array(atoms.get_tags()) #[1]*len(atoms) @@ -399,11 +399,11 @@ def write_lammps_data_water(self, atoms): for i in range(N_mol): fh.write('{:4d} {:4d} {:4d} {:4d}\n'.format(i*2+1,1,i*3+1,i*3+2)) fh.write('{:4d} {:4d} {:4d} {:4d}\n'.format(i*2+2,1,i*3+1,i*3+3)) - + fh.write('\nAngles \n\n') for i in range(N_angle): fh.write('{:4d} {:4d} {:4d} {:4d} {:4d}\n'.format(i+1,1,i*3+2,i*3+1,i*3+3)) - + def update(self, atoms): if not hasattr(self, 'atoms') or self.atoms != atoms: @@ -458,7 +458,7 @@ def __init__(self, temp=500): "unfix f1\n", "unfix f2\n", ] - + self.dict['GB_md'] = [ "# ---------- Run GB MD ---------\n", "thermo 100\n", diff --git a/pyxtal/representation.py b/pyxtal/representation.py index 5481c9e5..3c558692 100644 --- a/pyxtal/representation.py +++ b/pyxtal/representation.py @@ -243,7 +243,7 @@ def from_string(cls, inputs, smiles, composition=None): msg = "Composition is inconsistent: {:d}/{:d}\n".format(sum(composition), n_site) msg += str(inputs) raise ValueError(msg) - n_cell += 1 + #n_cell += 1 for i, smile in enumerate(smiles): if smile.endswith('.smi'): @@ -253,16 +253,20 @@ def from_string(cls, inputs, smiles, composition=None): n_mol = 4 else: n_torsion = len(find_rotor_from_smile(smile)) - n_mol = 7 + n_torsion + n_mol = 8 + n_torsion # (wp_id, x, y, z, ori_x, ori_y, ori_z, inv) + torsion #inversion - inputs[n_cell+n_mol-2] = int(inputs[n_cell+n_mol-2]) - x.append(inputs[n_cell-1:n_cell+n_mol-1]) + #print(n_mol, n_cell, len(inputs)) + inputs[n_cell] = int(inputs[n_cell]) + inputs[n_cell+n_mol-1] = int(inputs[n_cell+n_mol-1]) + x.append(inputs[n_cell:n_cell+n_mol])#; print('string', x[-1]) n_cell += n_mol + assert(n_cell == len(inputs)) return cls(x, smiles) def to_standard_setting(self): xtal = self.to_pyxtal() - rep0 = representation.from_pyxtal(xtal, standard=True) + xtal.optimize_lattice(standard=True) + rep0 = representation.from_pyxtal(xtal) self.x = rep0.x def to_pyxtal(self, smiles=None, composition=None): @@ -335,7 +339,7 @@ def to_pyxtal(self, smiles=None, composition=None): dicts['center'] = v[1:4] if smile not in ["Cl-"]: dicts['orientation'] = np.array(v[4:7]) - dicts['rotor'] = v[7:-1] + dicts['rotor'] = v[7:-1]#; print('ro', dicts['rotor']) dicts['reflect'] = int(v[-1]) site = mol_site.from_1D_dicts(dicts) @@ -389,7 +393,7 @@ def to_string(self, time=None, eng=None, tag=None): strs += "{:5.1f} ".format(c) # data for molecule - strs += "{:d} ".format(len(x)-1) + strs += "{:d} ".format(len(x)-1)#; print(x[1]) for i in range(1, len(x)): strs += "{:d} ".format(x[i][0]) for v in x[i][1:4]: @@ -472,7 +476,7 @@ def get_dist(self, rep): #aspirin smiles = ['CC(=O)OC1=CC=CC=C1C(=O)O'] - x = [[81,11.43,6.49,11.19,83.31],[0.77,0.57,0.53,48.55,24.31,145.94,-77.85,-4.40,170.86,False]] + x = [[81,11.43,6.49,11.19,83.31],[0, 0.77,0.57,0.53,48.55,24.31,145.94,-77.85,-4.40,170.86,False]] #rep0 = representation(x, smiles) #print(rep0.to_string()) rep1 = representation(x, smiles) @@ -481,16 +485,20 @@ def get_dist(self, rep): rep2 = representation.from_pyxtal(xtal) print(rep2.to_pyxtal()) print(rep2.to_string()) - string = "82 11.43 6.49 11.19 83.31 1 0.77 0.57 0.53 48.55 24.31 145.9 -77.85 -4.40 170.9 0" + + print('Test read from string') + string = "82 11.43 6.49 11.19 83.31 1 0 0.77 0.57 0.53 48.55 24.31 145.9 -77.85 -4.40 170.9 0" rep3 = representation.from_string(string, smiles) print(rep3.to_string()) print(rep3.to_pyxtal()) + #x = rep3.to_pyxtal(); x.optimize_lattice(standard=True); print(x) rep3.to_standard_setting() print(rep3.to_pyxtal()) print(rep3.to_string()) - string1 = "81 14.08 6.36 25.31 83.9 1 0.83 0.40 0.63 136.6 -21.6 -151.1 -101.1 -131.2 154.7 -176.4 -147.8 178.2 -179.1 -53.3 0" - string2 = "81 14.08 6.36 25.31 83.9 1 0.03 0.84 0.89 149.1 -8.0 -37.8 -39.9 -104.2 176.2 -179.6 137.8 -178.5 -173.3 -103.6 0" + print("Test other cases") + string1 = "81 14.08 6.36 25.31 83.9 1 0 0.83 0.40 0.63 136.6 -21.6 -151.1 -101.1 -131.2 154.7 -176.4 -147.8 178.2 -179.1 -53.3 0" + string2 = "81 14.08 6.36 25.31 83.9 1 0 0.03 0.84 0.89 149.1 -8.0 -37.8 -39.9 -104.2 176.2 -179.6 137.8 -178.5 -173.3 -103.6 0" smiles = ['CC1=CC=C(C=C1)S(=O)(=O)C2=C(N=C(S2)C3=CC=C(C=C3)NC(=O)OCC4=CC=CC=C4)C'] rep4 = representation.from_string(string1, smiles) rep5 = representation.from_string(string2, smiles)