diff --git a/qforce/data/default_md.mdp b/qforce/data/default_md.mdp new file mode 100644 index 0000000..b3f78e8 --- /dev/null +++ b/qforce/data/default_md.mdp @@ -0,0 +1,53 @@ +; MDP file + +; Run parameters +integrator = md ; leap-frog integrator +nsteps = 10000000 ; 10'000'000 = 3 ns +dt = 0.0001 ; in ps (0.001 = 1 fs) +nstcomm = 1 ; freq. for cm-motion removal +ld_seed = -1 + +; Bond constraints +continuation = no ; continue from npt equilibration +constraints = none ; constrain hydrogen bond lengths +constraint_algorithm = lincs ; default +lincs_order = 4 ; default + +; X/V/F/E outputs +nstxout = 50000 ; pos out --- 1000 ps +nstvout = 50000 ; vel out --- 1000 ps +nstfout = 0 ; force out --- no +nstlog = 50000 ; energies to log (20 ps) +nstenergy = 50000 ; energies to energy file +nstxout-compressed = 500 ; xtc, 1 ps +compressed-x-precision = 100000 + +; Neighbour list +cutoff-scheme = Verlet ; +ns_type = grid ; neighlist type +nstlist = 20 ; Freq. to update neighbour list +rlist = 1.2 ; nm (cutoff for short-range NL) + +; Coulomb interactions +coulombtype = Cut-off +rcoulomb = 1.1 ; nm (direct space sum cut-off) +;optimize_fft = yes ; optimal FFT plan for the grid + +; van der Waals interactions +vdwtype = Cut-off ; Van der Waals interactions +rvdw = 1.1 ; nm (LJ cut-off) +DispCorr = No ; use dispersion correction + +; Energy monitoring +energygrps = System + +; Temperature coupling is on +tcoupl = V-rescale ; modified Berendsen thermostat +tc-grps = MOL Argon ; two coupling groups - more accurate +tau_t = 1.0 1.0 ; time constant, in ps +ref_t = 273 273 ; reference temperature, one for each group, in K + +; Periodic boundary conditions +pbc = xyz ; 3-D PBC +; Velocity generation +gen_vel = yes \ No newline at end of file diff --git a/qforce/dihedral_scan.py b/qforce/dihedral_scan.py index 6f4a030..3db97b0 100644 --- a/qforce/dihedral_scan.py +++ b/qforce/dihedral_scan.py @@ -1,3 +1,4 @@ +# fmt: off import subprocess import os import shutil @@ -19,6 +20,8 @@ from .calculator import QForce from .forces import get_dihed, get_dist +from typing import List, Tuple + """ Fit all dihedrals togethers after the scans? @@ -46,7 +49,7 @@ class DihedralScan(Colt): break_co_bond = no :: bool # Method for doing the MM relaxed dihedral scan -method = qforce :: str :: [qforce, gromacs] +method = qforce :: str :: [qforce, gromacs, abf] # The executable for gromacs - necessary if scan method is gromacs gromacs_exec = gmx :: str @@ -71,6 +74,7 @@ class DihedralScan(Colt): """ def __init__(self, fragments, mol, job, all_config): + self.job = job self.frag_dir = job.frag_dir self.job_name = job.name self.mdp_file = f'{job.md_data}/default.mdp' @@ -79,6 +83,8 @@ def __init__(self, fragments, mol, job, all_config): self.scan = getattr(self, f'scan_dihed_{self.config.method.lower()}') self.move_capping_atoms(fragments) + print(f"---------- Scanning with {self.config.method.lower()} ----------") + fragments, all_dih_terms, weights = self.arrange_data(mol, fragments) final_energy, params = self.scan_dihedrals(fragments, mol, all_config, all_dih_terms, weights) @@ -163,18 +169,25 @@ def finalize_results(self, fragments, final_energy, all_dih_terms, params): print(' Please check manually to see if you find the accuracy satisfactory.\n') def scan_dihedrals(self, fragments, mol, all_config, all_dih_terms, weights): + removedir = [] + for n_run in range(self.config.n_dihed_scans): energy_diffs, md_energies = [], [] + for n_fit, frag in enumerate(fragments, start=1): print(f'Run {n_run+1}/{self.config.n_dihed_scans}, fitting dihedral ' f'{n_fit}/{len(fragments)}: {frag.id}') scan_dir = f'{self.frag_dir}/{frag.id}' - make_scan_dir(scan_dir) + + if frag.id not in removedir: + make_scan_dir(scan_dir) + removedir.append(frag.id) for term in frag.fit_terms: term['angles'] = [] + md_energy = self.scan(all_config, frag, scan_dir, mol, n_run) md_energy -= md_energy.min() @@ -262,7 +275,106 @@ def scan_dihed_gromacs(self, all_config, frag, scan_dir, mol, n_run): coords = read(f'{step_dir}/geom.gro').get_positions() self.calc_fit_angles(frag, coords) frag.coords[i] = coords + + # for i, term in enumerate(frag.fit_terms): + # print(f"Angle series (original, rad) {i}: {term['angles']}") + # print(f"Angle series (deg) {i}: {np.degrees(term['angles'])}") + return np.array(md_energies) + +# fmt: on + def scan_dihed_abf(self, all_config, frag, scan_dir, mol, n_run): + ff = ForceField( + self.job_name, + all_config, + frag, + frag.neighbors, + exclude_all=frag.remove_non_bonded, + ) + + # QM angles + step_size = all_config.qm.scan_step_size + first_qm_coord = frag.coords[0] + qm_term = frag.fit_terms[0] # term containing atom indices of the "main" scanned dihedral + first_qm_angle = np.degrees(get_dihed(first_qm_coord[qm_term['atomids']])[0]) + last_qm_angle = first_qm_angle + ((360/step_size)-1)*step_size + wraparound = first_qm_angle+180 + + # Managing MM angles + for term in frag.fit_terms: + angles = [] + for coord in frag.coords: + angle = get_dihed(coord[term['atomids']])[0] + angles.append(angle) + + term['angles'] = angles + + run_dir = f"{scan_dir}/abf_run_{n_run}" + make_scan_dir(run_dir) + + ff.write_gromacs(run_dir, frag, frag.coords[0], abf=True) + shutil.copy2(f'{self.job.md_data}/default_md.mdp', run_dir) + + index1, index2, index3, index4 = frag.scanned_atomids + + + colvar_string = f"""colvar {{ + name phi + width {step_size} + lowerBoundary {first_qm_angle} + upperBoundary {last_qm_angle} + # H C C H + dihedral {{ + wraparound {wraparound} + group1 {{ atomNumbers {index1+1}}} + group2 {{ atomNumbers {index2+1}}} + group3 {{ atomNumbers {index3+1}}} + group4 {{ atomNumbers {index4+1}}} + }} + extendedLagrangian on + extendedFluctuation 0.3 +}} +abf {{ + colvars phi + fullSamples 200 +}}""" + with open(f"{run_dir}/colvars.dat", "w") as colvarsfile: + colvarsfile.write(colvar_string) + + # create argon single atom pdb file on the fly + with open(f"{run_dir}/argon.pdb", "w") as argon_pdb_file: + argon_pdb_file.write("HETATM 1 Ar UNL 1 0.000 0.000 0.000 1.00 0.00 Ar") + + argon_itp_string = ''' +[ moleculetype ] +; Name nrexcl +Argon 3 + +[ atoms ] +; nr type resnr residue atom cgnr charge mass typeB chargeB massB +; residue 1 RES rtp RES q 0.0 + 1 opls_097 1 Argon Ar 1 0.00000000 39.948000 ; qtot 0.000000 +''' + with open(f"{run_dir}/argon.itp", "w") as argon_itp_file: + argon_itp_file.write(argon_itp_string) + + gromacs_exec = all_config.scan.gromacs_exec + + # create sample with gromacs insert + insert = subprocess.Popen([gromacs_exec, 'insert-molecules', '-f', 'gas.gro', '-ci', 'argon.pdb', '-nmol', '40', '-o', 'system.gro'], + cwd=run_dir, stdout=subprocess.PIPE, + stderr=subprocess.STDOUT) + insert.wait() + # run! + run_gromacs(run_dir, all_config.scan.gromacs_exec, ff.polar_title, abf = True) + + energy = np.genfromtxt(f"{run_dir}/em.pmf", comments='#', usecols=1) + + return energy + + + +# fmt: off @staticmethod def calc_fit_angles(frag, coords): @@ -383,6 +495,11 @@ def _set_symmetrize(self): 'direct': direct[i] == '+'}) return sym_dict +def wrap_angle(angle: float) -> float: + angle = angle % 360 + if angle > 180: + angle -= 360 + return angle def get_periodic_angle(angle): if angle > 360: @@ -393,6 +510,8 @@ def get_periodic_angle(angle): def get_periodic_angles(angles): + if isinstance(angles, list): + angles = np.array(angles) angles[angles > 360] -= 360 angles[angles < 0] += 360 return angles @@ -419,10 +538,18 @@ def make_scan_dir(scan_name): os.makedirs(scan_name) -def run_gromacs(directory, gromacs_exec, polar_title): +def run_gromacs(directory, gromacs_exec, polar_title, abf=False): + if abf: + mdp_file = 'default_md.mdp' + structure_file = 'system.gro' + colvars = ['-colvars', 'colvars.dat'] + else: + mdp_file = 'default.mdp' + structure_file = f'gas{polar_title}.gro' + colvars = [] attempt, returncode = 0, 1 - grompp = subprocess.Popen([gromacs_exec, 'grompp', '-f', 'default.mdp', '-p', - f'gas{polar_title}.top', '-c', f'gas{polar_title}.gro', '-o', + grompp = subprocess.Popen([gromacs_exec, 'grompp', '-f', mdp_file, '-p', + f'gas{polar_title}.top', '-c', structure_file, '-o', 'em.tpr', '-po', 'em.mdp', '-maxwarn', '10'], cwd=directory, stdout=subprocess.PIPE, stderr=subprocess.STDOUT) @@ -430,7 +557,7 @@ def run_gromacs(directory, gromacs_exec, polar_title): while returncode != 0 and attempt < 5: check_gromacs_termination(grompp) - mdrun = subprocess.Popen([gromacs_exec, 'mdrun', '-deffnm', 'em'], + mdrun = subprocess.Popen([gromacs_exec, 'mdrun', '-nt', '1', '-deffnm', 'em', *colvars], cwd=directory, stdout=subprocess.PIPE, stderr=subprocess.STDOUT) mdrun.wait() @@ -483,8 +610,10 @@ def calc_multi_rb_matrix(fragments, all_dih_terms, n_total_scans): matrix = np.zeros((n_total_scans, n_terms)) for frag in fragments: n_scans = len(frag.qm_angles) + for term in frag.fit_terms: term_idx = all_dih_terms.index(term['name']) + matrix[scan_sum:scan_sum+n_scans, term_idx*6:(term_idx*6)+6] += calc_rb(term['angles']) scan_sum += n_scans @@ -507,3 +636,4 @@ def calc_rb_pot(params, angles): for i in range(1, 6): rb += params[i]*np.cos(np.array(angles)-np.pi)**i return rb +# fmt: on diff --git a/qforce/forcefield.py b/qforce/forcefield.py index 2c71264..9bd1589 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -1,4 +1,6 @@ +# fmt: off import numpy as np + # from .elements import ATOM_SYM, ATOMMASS from .molecule.non_bonded import calc_sigma_epsilon @@ -6,7 +8,7 @@ from .misc import LOGO_SEMICOL -class ForceField(): +class ForceField: def __init__(self, job_name, config, mol, neighbors, exclude_all=[]): self.polar = config.ff._polar self.mol_name = job_name @@ -23,27 +25,42 @@ def __init__(self, job_name, config, mol, neighbors, exclude_all=[]): self.masses = [round(ATOMMASS[i], 5) for i in self.elements] self.exclusions = self.make_exclusions(mol.non_bonded, neighbors, exclude_all) self.pairs = self.make_pairs(neighbors, mol.non_bonded) + self.mol = mol + self.config = config if self.polar: - self.polar_title = '_polar' + self.polar_title = "_polar" else: - self.polar_title = '' - - def write_gromacs(self, directory, mol, coords): - self.write_itp(mol, directory) - self.write_top(directory) - self.write_gro(directory, coords, mol.non_bonded.alpha_map) + self.polar_title = "" + + def write_gromacs(self, directory, mol, coords, abf=False): + self.write_itp(mol, directory, abf) + self.write_top(directory, abf) + if abf: + box = [5.0, 5.0, 5.0] + self.write_gro(directory, coords, mol.non_bonded.alpha_map, box=box) + else: + self.write_gro(directory, coords, mol.non_bonded.alpha_map) - def write_top(self, directory): + def write_top(self, directory, abf=False): with open(f"{directory}/gas{self.polar_title}.top", "w") as top: # defaults top.write("\n[ defaults ]\n") top.write("; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ\n") - top.write(f"{1:>8} {self.comb_rule:>12} {'yes':>12} {self.fudge_lj:>12} " - f"{self.fudge_q:>12}\n\n\n") + top.write( + f"{1:>8} {self.comb_rule:>12} {'yes':>12} {self.fudge_lj:>12} " + f"{self.fudge_q:>12}\n\n\n" + ) top.write("; Include the molecule ITP\n") - top.write(f'#include "./{self.mol_name}_qforce{self.polar_title}.itp"\n\n\n') + top.write( + f'#include "./{self.mol_name}_qforce{self.polar_title}.itp"\n' + ) + + if abf: + top.write('#include "argon.itp"\n') + + top.write("\n\n") size = len(self.mol_name) top.write("[ system ]\n") @@ -53,16 +70,20 @@ def write_top(self, directory): top.write("[ molecules ]\n") top.write(f"; {' '*(size-10)}compound n_mol\n") top.write(f"{' '*(10-size)}{self.mol_name} 1\n") + if abf: + top.write("Argon 40\n") - def write_gro(self, directory, coords, alpha_map, box=[20., 20., 20.]): + def write_gro(self, directory, coords, alpha_map, box=[20.0, 20.0, 20.0]): n_atoms = self.n_atoms if self.polar: n_atoms += len(alpha_map.keys()) - coords_nm = coords*0.1 + coords_nm = coords * 0.1 with open(f"{directory}/gas{self.polar_title}.gro", "w") as gro: gro.write(f"{self.mol_name}\n") gro.write(f"{n_atoms:>6}\n") - for i, (a_name, coord) in enumerate(zip(self.atom_names, coords_nm), start=1): + for i, (a_name, coord) in enumerate( + zip(self.atom_names, coords_nm), start=1 + ): gro.write(f"{1:>5}{self.residue:<5}") gro.write(f"{a_name:>5}{i:>5}") gro.write(f"{coord[0]:>8.3f}{coord[1]:>8.3f}{coord[2]:>8.3f}\n") @@ -71,12 +92,14 @@ def write_gro(self, directory, coords, alpha_map, box=[20., 20., 20.]): gro.write(f"{2:>5}{self.residue:<5}{f'D{i}':>5}{drude+1:>5}") gro.write(f"{coords_nm[atom][0]:>8.3f}{coords_nm[atom][1]:>8.3f}") gro.write(f"{coords_nm[atom][2]:>8.3f}\n") - gro.write(f'{box[0]:>12.5f}{box[1]:>12.5f}{box[2]:>12.5f}\n') + gro.write(f"{box[0]:>12.5f}{box[1]:>12.5f}{box[2]:>12.5f}\n") - def write_itp(self, mol, directory): - with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.itp", "w") as itp: + def write_itp(self, mol, directory, abf): + with open( + f"{directory}/{self.mol_name}_qforce{self.polar_title}.itp", "w" + ) as itp: itp.write(LOGO_SEMICOL) - self.write_itp_atoms_and_molecule(itp, mol.non_bonded) + self.write_itp_atoms_and_molecule(itp, mol.non_bonded, abf) if self.polar: self.write_itp_polarization(itp, mol.non_bonded) self.write_itp_bonds(itp, mol.terms, mol.non_bonded.alpha_map) @@ -84,7 +107,7 @@ def write_itp(self, mol, directory): self.write_itp_dihedrals(itp, mol.terms) self.write_itp_pairs(itp) self.write_itp_exclusions(itp) - itp.write('\n') + itp.write("\n") def convert_to_gromacs_nonbonded(self, non_bonded): a_types, nb_pairs, nb_1_4 = {}, {}, {} @@ -117,8 +140,10 @@ def convert_to_gromacs_nonbonded(self, non_bonded): return a_types, nb_pairs, nb_1_4 - def write_itp_atoms_and_molecule(self, itp, non_bonded): - gro_atomtypes, gro_nonbonded, gro_1_4 = self.convert_to_gromacs_nonbonded(non_bonded) + def write_itp_atoms_and_molecule(self, itp, non_bonded, abf): + gro_atomtypes, gro_nonbonded, gro_1_4 = self.convert_to_gromacs_nonbonded( + non_bonded + ) # atom types itp.write("\n[ atomtypes ]\n") @@ -129,7 +154,10 @@ def write_itp_atoms_and_molecule(self, itp, non_bonded): for lj_type, lj_params in gro_atomtypes.items(): itp.write(f'{lj_type:>8} {non_bonded.lj_atomic_number[lj_type]:>7} {0:>8.4f} {0:>8.4f} {"A":>5} ') - itp.write(f'{lj_params[0]:>12.5e} {lj_params[1]:>12.5e}\n') + itp.write(f"{lj_params[0]:>12.5e} {lj_params[1]:>12.5e}\n") + + if abf: + itp.write("opls_097 18 0.0000 0.0000 A 3.40100e-01 9.78638e-01") if self.polar: itp.write(f'{"DP":>8} {0:>8.4f} {0:>8.4f} {"S":>2} {0:>12.5e} {0:>12.5e}\n') @@ -142,8 +170,8 @@ def write_itp_atoms_and_molecule(self, itp, non_bonded): itp.write("; name1 name2 type sigma epsilon\n") for lj_types, lj_params in gro_nonbonded.items(): - itp.write(f'{lj_types[0]:>8} {lj_types[1]:>8} 1') - itp.write(f'{lj_params[0]:>15.5e} {lj_params[1]:>15.5e}\n') + itp.write(f"{lj_types[0]:>8} {lj_types[1]:>8} 1") + itp.write(f"{lj_params[0]:>15.5e} {lj_params[1]:>15.5e}\n") # Write 1-4 pair types if self.n_excl == 2 and gro_1_4 != {}: @@ -153,11 +181,11 @@ def write_itp_atoms_and_molecule(self, itp, non_bonded): else: itp.write("; name1 name2 type sigma epsilon\n") for lj_types, lj_params in gro_1_4.items(): - itp.write(f'{lj_types[0]:>8} {lj_types[1]:>8} 1') - itp.write(f'{lj_params[0]:>15.5e} {lj_params[1]:>15.5e}\n') + itp.write(f"{lj_types[0]:>8} {lj_types[1]:>8} 1") + itp.write(f"{lj_params[0]:>15.5e} {lj_params[1]:>15.5e}\n") # molecule type - space = " "*(len(self.mol_name)-5) + space = " " * (len(self.mol_name) - 5) itp.write("\n[ moleculetype ]\n") itp.write(f";{space}name nrexcl\n") itp.write(f"{self.mol_name}{3:>7}\n") @@ -165,22 +193,23 @@ def write_itp_atoms_and_molecule(self, itp, non_bonded): # atoms itp.write("\n[ atoms ]\n") itp.write("; nr type resnr resnm atom cgrp charge mass\n") - for i, (lj_type, a_name, q, mass) in enumerate(zip(non_bonded.lj_types, self.atom_names, - self.q, self.masses), start=1): - itp.write(f'{i:>5} {lj_type:>9} {1:>6} {self.residue:>6} {a_name:>7} {i:>5} ') - itp.write(f'{q:>11.5f} {mass:>10.5f}\n') + for i, (lj_type, a_name, q, mass) in enumerate( + zip(non_bonded.lj_types, self.atom_names, self.q, self.masses), start=1 + ): + itp.write(f"{i:>5} {lj_type:>9} {1:>6} {self.residue:>6} {a_name:>7} {i:>5} ") + itp.write(f"{q:>11.5f} {mass:>10.5f}\n") if self.polar: for i, (atom, drude) in enumerate(non_bonded.alpha_map.items(), start=1): itp.write(f'{drude+1:>5} {"DP":>9} {2:>6} {"DRU":>6} {f"D{i}":>7} {atom+1:>5}') - itp.write(f'{-8.:>11.5f} {0.:>10.5f}\n') + itp.write(f"{-8.:>11.5f} {0.:>10.5f}\n") def write_itp_polarization(self, itp, non_bonded): # polarization itp.write("\n[ polarization ]\n") itp.write("; i j f alpha\n") for atom, drude in non_bonded.alpha_map.items(): - alpha = non_bonded.alpha[atom]*1e-3 + alpha = non_bonded.alpha[atom] * 1e-3 itp.write(f"{atom+1:>6} {drude+1:>6} {1:>6} {alpha:>14.8f}\n") # # thole polarization @@ -201,99 +230,117 @@ def write_itp_pairs(self, itp): def write_itp_bonds(self, itp, terms, alpha_map): itp.write("\n[ bonds ]\n") itp.write("; ai aj f r0 k_bond\n") - for bond in terms['bond']: + for bond in terms["bond"]: ids = bond.atomids + 1 equ = bond.equ * 0.1 fconst = bond.fconst * 100 - itp.write(f'{ids[0]:>6} {ids[1]:>6} {1:>6} {equ:>10.5f} {fconst:>10.0f}\n') + itp.write(f"{ids[0]:>6} {ids[1]:>6} {1:>6} {equ:>10.5f} {fconst:>10.0f}\n") if self.polar: - itp.write('; ai aj f - polar connections\n') - for bond in terms['bond']: + itp.write("; ai aj f - polar connections\n") + for bond in terms["bond"]: a1, a2 = bond.atomids if a2 in alpha_map.keys(): - itp.write(f'{a1+1:>6} {alpha_map[a2]+1:>6} {5:>6}\n') + itp.write(f"{a1+1:>6} {alpha_map[a2]+1:>6} {5:>6}\n") if a1 in alpha_map.keys(): - itp.write(f'{a2+1:>6} {alpha_map[a1]+1:>6} {5:>6}\n') + itp.write(f"{a2+1:>6} {alpha_map[a1]+1:>6} {5:>6}\n") def write_itp_angles(self, itp, terms): itp.write("\n[ angles ]\n") itp.write("; ai aj ak f theta0 k_theta") if self.urey: - itp.write(' r0 k_bond') - itp.write('\n') + itp.write(" r0 k_bond") + itp.write("\n") - for angle in terms['angle']: + for angle in terms["angle"]: ids = angle.atomids + 1 equ = np.degrees(angle.equ) fconst = angle.fconst if self.urey: - urey = [term for term in terms['urey'] if np.array_equal(term.atomids, - angle.atomids)] + urey = [ + term + for term in terms["urey"] + if np.array_equal(term.atomids, angle.atomids) + ] if not self.urey or len(urey) == 0: - itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {1:>6} {equ:>10.3f} ' - f'{fconst:>10.3f}\n') + itp.write(f"{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {1:>6} {equ:>10.3f} " + f"{fconst:>10.3f}\n") else: urey_equ = urey[0].equ * 0.1 urey_fconst = urey[0].fconst * 100 - itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {5:>6} {equ:>10.3f} ' - f'{fconst:>10.3f} {urey_equ:>10.5f} {urey_fconst:>10.1f}\n') + itp.write( + f"{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {5:>6} {equ:>10.3f} " + f"{fconst:>10.3f} {urey_equ:>10.5f} {urey_fconst:>10.1f}\n" + ) def write_itp_dihedrals(self, itp, terms): - if len(terms['dihedral']) > 0: + if len(terms["dihedral"]) > 0: itp.write("\n[ dihedrals ]\n") # rigid dihedrals - if len(terms['dihedral/rigid']) > 0: + if len(terms["dihedral/rigid"]) > 0: itp.write("; rigid dihedrals \n") itp.write("; ai aj ak al f theta0 k_theta\n") - for dihed in terms['dihedral/rigid']: + for dihed in terms["dihedral/rigid"]: ids = dihed.atomids + 1 equ = np.degrees(dihed.equ) fconst = dihed.fconst - itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {2:>6} {equ:>11.3f} ') - itp.write(f'{fconst:>13.3f}\n') + itp.write(f"{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {2:>6} {equ:>11.3f} ") + itp.write(f"{fconst:>13.3f}\n") # improper dihedrals - if len(terms['dihedral/improper']) > 0: + if len(terms["dihedral/improper"]) > 0: itp.write("; improper dihedrals \n") itp.write("; ai aj ak al f theta0 k_theta\n") - for dihed in terms['dihedral/improper']: + for dihed in terms["dihedral/improper"]: ids = dihed.atomids + 1 equ = np.degrees(dihed.equ) fconst = dihed.fconst - itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {2:>6} {equ:>11.3f} ') - itp.write(f'{fconst:>13.3f}\n') + itp.write( + f"{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {2:>6} {equ:>11.3f} " + ) + itp.write(f"{fconst:>13.3f}\n") # flexible dihedrals - if len(terms['dihedral/flexible']) > 0: + if len(terms["dihedral/flexible"]) > 0: itp.write("; flexible dihedrals \n") - itp.write('; ai aj ak al f c0 c1 c2') - itp.write(' c3 c4 c5\n') - - for dihed in terms['dihedral/flexible']: + itp.write("; ai aj ak al f c0 c1 c2") + itp.write(" c3 c4 c5\n") + + for dihed in terms["dihedral/flexible"]: + # print("\033[31m From forcefield.py:ForceField:write_itp_dihedral: \033[0m") + # print(type(dihed)) + # print(self.config.terms) + # print(self.mol.name, self.mol.dir) + # print(dihed.name, dihed.idx, dihed.atomids) + # print("Type of dihed.equ:") + # print(type(dihed.equ)) + # print("Contents of dihed.equ:") + # print(dihed.equ) ids = dihed.atomids + 1 c = dihed.equ - itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {3:>6} {c[0]:>11.3f} ') - itp.write(f'{c[1]:>11.3f} {c[2]:>11.3f} {c[3]:>11.3f} {c[4]:>11.3f} {c[5]:>11.3f}\n') + itp.write(f"{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {3:>6} {c[0]:>11.3f} ") + itp.write(f"{c[1]:>11.3f} {c[2]:>11.3f} {c[3]:>11.3f} {c[4]:>11.3f} {c[5]:>11.3f}\n") # inversion dihedrals - if len(terms['dihedral/inversion']) > 0: + if len(terms["dihedral/inversion"]) > 0: itp.write("; inversion dihedrals \n") - itp.write('; ai aj ak al f c0 c1 c2\n') + itp.write("; ai aj ak al f c0 c1 c2\n") - for dihed in terms['dihedral/inversion']: + for dihed in terms["dihedral/inversion"]: ids = dihed.atomids + 1 c0, c1, c2 = convert_to_inversion_rb(dihed.fconst, dihed.equ) - itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {3:>6}' - f'{c0:>11.3f} {c1:>11.3f} {c2:>11.3f} {0:>11.1f} {0:>11.1f} {0:>11.1f}\n') + itp.write( + f"{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {3:>6}" + f"{c0:>11.3f} {c1:>11.3f} {c2:>11.3f} {0:>11.1f} {0:>11.1f} {0:>11.1f}\n" + ) def write_itp_exclusions(self, itp): if any(len(exclusion) > 0 for exclusion in self.exclusions): @@ -301,8 +348,8 @@ def write_itp_exclusions(self, itp): for i, exclusion in enumerate(self.exclusions): if len(exclusion) > 0: exclusion = sorted(set(exclusion)) - itp.write("{} ".format(i+1)) - itp.write(("{} "*len(exclusion)).format(*exclusion)) + itp.write("{} ".format(i + 1)) + itp.write(("{} " * len(exclusion)).format(*exclusion)) itp.write("\n") def make_pairs(self, neighbors, non_bonded): @@ -315,47 +362,54 @@ def make_pairs(self, neighbors, non_bonded): polar_pairs.append([a1, non_bonded.alpha_map[a2]]) if a1 in non_bonded.alpha_map.keys(): polar_pairs.append([a2, non_bonded.alpha_map[a1]]) - if a1 in non_bonded.alpha_map.keys() and a2 in non_bonded.alpha_map.keys(): + if (a1 in non_bonded.alpha_map.keys() and a2 in non_bonded.alpha_map.keys()): polar_pairs.append([non_bonded.alpha_map[a1], non_bonded.alpha_map[a2]]) - return non_bonded.pairs+polar_pairs + return non_bonded.pairs + polar_pairs def make_exclusions(self, non_bonded, neighbors, exclude_all): exclusions = [[] for _ in range(self.n_atoms)] # input exclusions for exclusions if outside of n_excl - for a1, a2 in non_bonded.exclusions+non_bonded.pairs: - if all([a2 not in neighbors[i][a1] for i in range(self.n_excl+1)]): - exclusions[a1].append(a2+1) + for a1, a2 in non_bonded.exclusions + non_bonded.pairs: + if all([a2 not in neighbors[i][a1] for i in range(self.n_excl + 1)]): + exclusions[a1].append(a2 + 1) # fragment capping atom exclusions for i in exclude_all: - exclusions[i].extend(np.arange(1, self.n_atoms+1)) + exclusions[i].extend(np.arange(1, self.n_atoms + 1)) if self.polar: - exclusions = self.polarize_exclusions(non_bonded.alpha_map, non_bonded.exclusions, - neighbors, exclude_all, exclusions) + exclusions = self.polarize_exclusions( + non_bonded.alpha_map, + non_bonded.exclusions, + neighbors, + exclude_all, + exclusions, + ) return exclusions - def polarize_exclusions(self, alpha_map, input_exclusions, neighbors, exclude_all, exclusions): + def polarize_exclusions( + self, alpha_map, input_exclusions, neighbors, exclude_all, exclusions + ): n_polar_atoms = len(alpha_map.keys()) exclusions.extend([[] for _ in range(n_polar_atoms)]) # input exclusions for exclu in input_exclusions: if exclu[0] in alpha_map.keys(): - exclusions[alpha_map[exclu[0]]].append(exclu[1]+1) + exclusions[alpha_map[exclu[0]]].append(exclu[1] + 1) if exclu[1] in alpha_map.keys(): - exclusions[alpha_map[exclu[1]]].append(exclu[0]+1) + exclusions[alpha_map[exclu[1]]].append(exclu[0] + 1) if exclu[0] in alpha_map.keys() and exclu[1] in alpha_map.keys(): - exclusions[alpha_map[exclu[0]]].append(alpha_map[exclu[1]]+1) + exclusions[alpha_map[exclu[0]]].append(alpha_map[exclu[1]] + 1) # fragment capping atom exclusions for i in exclude_all: - exclusions[i].extend(np.arange(self.n_atoms+1, self.n_atoms+n_polar_atoms+1)) + exclusions[i].extend(np.arange(self.n_atoms + 1, self.n_atoms + n_polar_atoms + 1)) if i in alpha_map.keys(): - exclusions[alpha_map[i]].extend(np.arange(1, self.n_atoms+n_polar_atoms+1)) + exclusions[alpha_map[i]].extend(np.arange(1, self.n_atoms + n_polar_atoms + 1)) return exclusions @@ -369,7 +423,7 @@ def get_atom_names(self): atom_dict[sym] = 1 else: atom_dict[sym] += 1 - atom_names.append(f'{sym}{atom_dict[sym]}') + atom_names.append(f"{sym}{atom_dict[sym]}") return atom_names def add_restraints(self, restraints, directory, fc=1000): @@ -377,9 +431,9 @@ def add_restraints(self, restraints, directory, fc=1000): itp.write("[ dihedral_restraints ]\n") itp.write("; ai aj ak al type phi dp kfac\n") for restraint in restraints: - a1, a2, a3, a4 = restraint[0]+1 + a1, a2, a3, a4 = restraint[0] + 1 phi = np.degrees(restraint[1]) - itp.write(f'{a1:>5} {a2:>5} {a3:>5} {a4:>5} {1:>5} {phi:>10.4f} 0.0 {fc}\n') + itp.write(f"{a1:>5} {a2:>5} {a3:>5} {a4:>5} {1:>5} {phi:>10.4f} 0.0 {fc}\n") def set_charge(self, non_bonded): q = np.copy(non_bonded.q) @@ -473,3 +527,4 @@ def set_charge(self, non_bonded): # elif in_section == "pairs": # self.pairs.append(sorted([int(line[0]), int(line[1])])) + # fmt: on diff --git a/qforce/qm/orca.py b/qforce/qm/orca.py index 37152dc..c5efd7d 100644 --- a/qforce/qm/orca.py +++ b/qforce/qm/orca.py @@ -5,6 +5,7 @@ import numpy as np from ase.units import Hartree, mol, kJ, Bohr from ase.io import read + # from .qm_base import WriteABC, ReadABC from ..elements import ATOM_SYM @@ -34,20 +35,22 @@ class Orca(Colt): """ - _method = ['qm_method_hessian', 'qm_method_opt', 'qm_method_charge', 'qm_method_sp'] + _method = ["qm_method_hessian", "qm_method_opt", "qm_method_charge", "qm_method_sp"] def __init__(self): - self.required_hessian_files = {'out_file': ['.out', '.log'], - 'hess_file': ['_opt.hess'], - 'pc_file': ['_charge.pc_chelpg'], - 'coord_file': ['_opt.xyz']} + self.required_hessian_files = { + "out_file": [".out", ".log"], + "hess_file": ["_opt.hess"], + "pc_file": ["_charge.pc_chelpg"], + "coord_file": ["_opt.xyz"], + } self.read = ReadORCA self.write = WriteORCA class WriteORCA(WriteABC): def hessian(self, file, job_name, config, coords, atnums): - """ Write the input file for hessian and charge calculation. + """Write the input file for hessian and charge calculation. Parameters ---------- @@ -66,45 +69,43 @@ def hessian(self, file, job_name, config, coords, atnums): # Write the coordinates file.write(f"* xyz {config.charge} {config.multiplicity}\n") self._write_coords(atnums, coords, file) - file.write(' *\n\n') + file.write(" *\n\n") file.write(f"%pal nprocs {config.n_proc} end\n") # ORCA uses MPI parallelization and a factor of 0.75 is used to # avoid ORCA using more than it is available. - file.write('%maxcore {}\n\n'.format(int( - config.memory / config.n_proc * 0.75))) + file.write("%maxcore {}\n\n".format(int(config.memory / config.n_proc * 0.75))) # Start compound job - file.write('%Compound\n\n') + file.write("%Compound\n\n") # Do the initial optimisation - file.write('New_Step\n') + file.write("New_Step\n") file.write(f"! opt {config.qm_method_opt} nopop\n") file.write(f'%base "{job_name}_initial"\n') - file.write('STEP_END\n\n') + file.write("STEP_END\n\n") # Do the hessian calculation - file.write('New_Step\n') + file.write("New_Step\n") file.write(f"! opt freq {config.qm_method_hessian} PModel nopop\n") file.write(f'%base "{job_name}_opt"\n') - file.write('STEP_END\n\n') + file.write("STEP_END\n\n") # Do the nbo calculation - file.write('New_Step\n') + file.write("New_Step\n") file.write(f"! {config.qm_method_hessian}\n") file.write(f'%base "{job_name}_nbo"\n') - file.write('%method\n MAYER_BONDORDERTHRESH 0\nend\n') - file.write('STEP_END\n\n') + file.write("%method\n MAYER_BONDORDERTHRESH 0\nend\n") + file.write("STEP_END\n\n") # Write the charge calculation input - file.write('New_Step\n') + file.write("New_Step\n") file.write(f"! {config.qm_method_charge} chelpg Hirshfeld nopop\n") file.write(f'%base "{job_name}_charge"\n') - file.write('STEP_END\n\n') - file.write('END\n') + file.write("STEP_END\n\n") + file.write("END\n") - def scan(self, file, job_name, config, coords, atnums, scanned_atoms, start_angle, charge, - multiplicity): - """ Write the input file for the dihedral scan and charge calculation. + def scan(self, file, job_name, config, coords, atnums, scanned_atoms, start_angle, charge, multiplicity): + """Write the input file for the dihedral scan and charge calculation. Parameters ---------- @@ -131,56 +132,53 @@ def scan(self, file, job_name, config, coords, atnums, scanned_atoms, start_angl # Write the coordinates file.write(f"* xyz {charge} {multiplicity}\n") self._write_coords(atnums, coords, file) - file.write(' *\n') + file.write(" *\n") file.write(f"%pal nprocs {config.n_proc} end\n") # ORCA uses MPI parallelization and a factor of 0.75 is used to # avoid ORCA using more than it is available. - file.write('%maxcore {}\n\n'.format(int( - config.memory / config.n_proc * 0.75))) + file.write("%maxcore {}\n\n".format(int(config.memory / config.n_proc * 0.75))) # Start compound job - file.write('%Compound\n\n') + file.write("%Compound\n\n") # Do the initial optimisation - file.write('New_Step\n') + file.write("New_Step\n") file.write(f"! opt {config.qm_method_opt} nopop\n") file.write(f'%base "{job_name}_opt"\n') self._write_constrained_atoms(file, scanned_atoms) - file.write('STEP_END\n\n') + file.write("STEP_END\n\n") # Get charge first - file.write('New_Step\n') + file.write("New_Step\n") # PModel used for initial guess such that using XTB would not pose a # problem. file.write(f"! {config.qm_method_charge} chelpg Hirshfeld PModel nopop\n") file.write(f'%base "{job_name}_charge"\n') - file.write('STEP_END\n\n') + file.write("STEP_END\n\n") # Do the scan - file.write('New_Step\n') + file.write("New_Step\n") file.write(f"! opt {config.qm_method_opt} nopop\n") file.write(f'%base "{job_name}_scan"\n') self._write_scanned_atoms(file, scanned_atoms, start_angle, config.scan_step_size) file.write(f"*xyzfile {charge} {multiplicity} {job_name}_opt.xyz\n") - file.write('STEP_END\n\n') + file.write("STEP_END\n\n") # Do the single point energy - file.write('New_Step\n') + file.write("New_Step\n") # PModel used for initial guess such that using XTB would not pose a # problem. file.write(f"! {config.qm_method_sp} PModel nopop\n") file.write(f'%base "{job_name}_sp"\n') - file.write( - f"*xyzfile {charge} {multiplicity} " - f"{job_name}_scan.allxyz\n") - file.write('STEP_END\n') + file.write(f"*xyzfile {charge} {multiplicity} " f"{job_name}_scan.allxyz\n") + file.write("STEP_END\n") # Close compound block - file.write('END\n') + file.write("END\n") @staticmethod def _write_scanned_atoms(file, scanned_atoms, start_angle, step_size): - """ Write the input line for dihedral scan. + """Write the input line for dihedral scan. Parameters ---------- @@ -200,15 +198,13 @@ def _write_scanned_atoms(file, scanned_atoms, start_angle, step_size): end_angle = start_angle + 360 - step_size n_steps = int(np.ceil(360 / step_size)) file.write("%geom Scan\n") - file.write(f"D {a1} {a2} {a3} {a4} = {start_angle:.2f}," - f" {end_angle:.2f}," - f" {n_steps}\n") + file.write(f"D {a1} {a2} {a3} {a4} = {start_angle:.2f}," f" {end_angle:.2f}," f" {n_steps}\n") file.write("end\n") file.write("end\n") @staticmethod def _write_constrained_atoms(file, scanned_atoms): - """ Write the input line for the constrained dihedral. + """Write the input line for the constrained dihedral. Parameters ---------- @@ -227,7 +223,7 @@ def _write_constrained_atoms(file, scanned_atoms): @staticmethod def _write_coords(atnums, coords, file): - """ Write the input coordinates. + """Write the input coordinates. Parameters ---------- @@ -240,13 +236,13 @@ def _write_coords(atnums, coords, file): """ for atnum, coord in zip(atnums, coords): elem = ATOM_SYM[atnum] - file.write(f'{elem :>3s} {coord[0]:>12.6f} {coord[1]:>12.6f} {coord[2]:>12.6f}\n') + file.write(f"{elem :>3s} {coord[0]:>12.6f} {coord[1]:>12.6f} {coord[2]:>12.6f}\n") class ReadORCA(ReadABC): @staticmethod def _read_orca_hess(hess_file): - """ Read the hessian matrix. + """Read the hessian matrix. For ORCA jobs, the file contain the hessian information is stored with the extension of hess. @@ -262,19 +258,19 @@ def _read_orca_hess(hess_file): An array of float of the size of ((n_atoms*3)**2+n_atoms*3)/2), which is the lower triangle of the hessian matrix. Unit： kJ/mol """ - with open(hess_file, 'r') as f: + with open(hess_file, "r") as f: text = f.read() - text = text[text.index('$hessian'):] - lines = text.split('\n')[1:] + text = text[text.index("$hessian") :] + lines = text.split("\n")[1:] n_atoms_times_3 = int(lines[0]) lines = lines[1:] hessian = np.empty((n_atoms_times_3, n_atoms_times_3)) for _ in range(int(np.ceil(n_atoms_times_3 / 5))): - trunk = lines[:(n_atoms_times_3 + 1)] - lines = lines[(n_atoms_times_3 + 1):] + trunk = lines[: (n_atoms_times_3 + 1)] + lines = lines[(n_atoms_times_3 + 1) :] cols = [int(col) for col in trunk[0].split()] for row in range(n_atoms_times_3): row_idx, *points = trunk[1 + row].split() @@ -293,7 +289,7 @@ def _read_orca_hess(hess_file): @staticmethod def _read_orca_esp(pc_file): - """ Read the point charge file. + """Read the point charge file. For ORCA jobs, the point charge computed by chelpg method will yield a file contain the charge information stored with the extension @@ -311,17 +307,17 @@ def _read_orca_esp(pc_file): point_charges : float A list of float of the size of n_atoms. """ - with open(pc_file, 'r') as f: - lines = f.read().split('\n') + with open(pc_file, "r") as f: + lines = f.read().split("\n") n_atoms = int(lines[0]) point_charges = [] - for line in lines[2: 2+n_atoms]: + for line in lines[2 : 2 + n_atoms]: point_charges.append(float(line.split()[1])) return n_atoms, point_charges @staticmethod def _read_orca_cm5(out_file): - """ Read the HIRSHFELD charge file. + """Read the HIRSHFELD charge file. Parameters ---------- @@ -335,18 +331,18 @@ def _read_orca_cm5(out_file): point_charges : float A list of float of the size of n_atoms. """ - file = open(out_file, 'r') + file = open(out_file, "r") line = file.readline() # Skip to HIRSHFELD ANALYSIS - while 'HIRSHFELD ANALYSIS' not in line: + while "HIRSHFELD ANALYSIS" not in line: line = file.readline() - while 'ATOM CHARGE SPIN' not in line: + while "ATOM CHARGE SPIN" not in line: line = file.readline() charges = [] - while 'TOTAL' not in line: + while "TOTAL" not in line: line = file.readline() if len(line.split()) == 4: atom_id, element, charge, _ = line.split() @@ -356,7 +352,7 @@ def _read_orca_cm5(out_file): @staticmethod def _read_orca_coords(coord_text): - """ Read the optimised coordinate string. + """Read the optimised coordinate string. The string is in the format of the standard xyz file. @@ -374,11 +370,11 @@ def _read_orca_coords(coord_text): coords : array An array of float of the shape (n_atoms, 3). """ - lines = coord_text.split('\n') + lines = coord_text.split("\n") n_atoms = int(lines[0]) coords = np.empty((n_atoms, 3)) elements = [] - for i, line in enumerate(lines[2: 2+n_atoms]): + for i, line in enumerate(lines[2 : 2 + n_atoms]): element, x, y, z = line.split() elements.append(ATOM_SYM.index(element)) coords[i, :] = (x, y, z) @@ -386,7 +382,7 @@ def _read_orca_coords(coord_text): @staticmethod def _read_orca_xyz(coord_file): - """ Read the optimised coordinate xyz file. + """Read the optimised coordinate xyz file. For ORCA jobs, the optimised geometry will be stored as a file with the extension xyz. @@ -413,7 +409,7 @@ def _read_orca_xyz(coord_file): @staticmethod def _read_orca_allxyz(coord_file): - """ Read the optimised coordinate allxyz file. + """Read the optimised coordinate allxyz file. For ORCA jobs, the scan optimised geometry will be stored as a file with the extension allxyz. @@ -433,17 +429,41 @@ def _read_orca_allxyz(coord_file): A list of array of float. The list has the length of the number of steps and the array has the shape of (n_atoms, 3). """ - with open(coord_file, 'r') as f: + # with open(coord_file, "r") as f: + # coord_text = f.read() + # coords = [] + # for text in coord_text.split(">\n"): + # n_atoms, elements, coord = ReadORCA._read_orca_coords(text) + # coords.append(coord) + # return n_atoms, elements, coords + + # Updated code to allow reading general xyz trajectories (structures not separated by > like in ORCA) + with open(coord_file, "r") as f: coord_text = f.read() coords = [] - for text in coord_text.split('>\n'): - n_atoms, elements, coord = ReadORCA._read_orca_coords(text) - coords.append(coord) + + block = [] + + counter = 0 + for line in coord_text.splitlines(): + if line.strip("\n").strip().isdigit() and block: + counter += 1 + n_atoms, elements, coord = ReadORCA._read_orca_coords("\n".join(block)) + coords.append(coord) + block = [line] + continue + if ">" in line: + continue + block.append(line) + + n_atoms, elements, coord = ReadORCA._read_orca_coords("\n".join(block)) + coords.append(coord) + return n_atoms, elements, coords @staticmethod def _read_orca_bond_order(out_file, n_atoms): - """ Read the bond order analysis. + """Read the bond order analysis. Parameters ---------- @@ -458,21 +478,22 @@ def _read_orca_bond_order(out_file, n_atoms): A list (length: n_atoms) of list (length: n_atoms) of float. representing the bond order between each atom pair. """ - item_match = re.compile('^\(\s*(\d+)-\w{1,2}\s*,\s*(\d+)-\w{1,2}\s*\)\s*:\s*(-?\w.+)$') - b_orders = [[0, ] * n_atoms for _ in range(n_atoms)] + # Matches patterns like B( 0-Cl, 1-C ) : 1.0623 + item_match = re.compile("^\(\s*(\d+)-\w{1,2}\s*,\s*(\d+)-\w{1,2}\s*\)\s*:\s*(-?\w.+)$") + b_orders = [[0] * n_atoms for _ in range(n_atoms)] - with open(out_file, 'r') as file: + with open(out_file, "r") as file: line = file.readline() # Skip to the step after geometry optimisation - while 'Mayer bond orders larger than' not in line: + while "Mayer bond orders larger than" not in line: line = file.readline() - + line = file.readline() while "-------" not in line: - items = line.split('B') + items = line.split("B(") for item in items: if item.strip(): - _m = re.match(item_match, item) + _m = re.match(item_match, "(" + item) i = int(_m.group(1)) j = int(_m.group(2)) bond_order = float(_m.group(3)) @@ -484,7 +505,7 @@ def _read_orca_bond_order(out_file, n_atoms): @staticmethod def _read_orca_dat(dat): - """ Read the scan parameter and energy from ORCA relaxscanact.dat or + """Read the scan parameter and energy from ORCA relaxscanact.dat or xyzact.dat file. Parameters @@ -499,17 +520,17 @@ def _read_orca_dat(dat): energies : list A list (length: steps) of the energy. """ - with open(dat, 'r') as f: + with open(dat, "r") as f: coord_text = f.read().strip() parameter = [] energies = [] - for line in coord_text.split('\n'): + for line in coord_text.split("\n"): parameter.append(float(line.split()[0])) energies.append(float(line.split()[1])) return parameter, energies def hessian(self, config, out_file, hess_file, pc_file, coord_file): - """ Extract information from all the relevant files. + """Extract information from all the relevant files. Parameters ---------- @@ -554,10 +575,10 @@ def hessian(self, config, out_file, hess_file, pc_file, coord_file): charge = config.charge multiplicity = config.multiplicity b_orders = self._read_orca_bond_order(out_file, n_atoms) - return n_atoms, charge, multiplicity, elements, coords, hessian, b_orders, point_charges + return (n_atoms, charge, multiplicity, elements, coords, hessian, b_orders, point_charges) def scan(self, config, file_name): - """ Read data from the scan file. + """Read data from the scan file. Parameters ---------- @@ -583,17 +604,19 @@ def scan(self, config, file_name): """ base, ext = os.path.splitext(file_name) point_charges = {} + if config.charge_method == "cm5": n_atoms, charges = self._read_orca_cm5(file_name) + point_charges["cm5"] = charges + if config.charge_method == "esp": - n_atoms, charges = self._read_orca_esp( - '{}_charge.pc_chelpg'.format(base)) + n_atoms, charges = self._read_orca_esp("{}_charge.pc_chelpg".format(base)) point_charges["esp"] = charges - n_atoms, elements, coords = self._read_orca_allxyz( - '{}_scan.allxyz'.format(base)) - angles, energies = self._read_orca_dat('{}_scan.relaxscanact.dat'.format(base)) - if os.path.isfile('{}_sp.xyzact.dat'.format(base)): - _, energies = self._read_orca_dat('{}_sp.xyzact.dat'.format(base)) + n_atoms, elements, coords = self._read_orca_allxyz("{}_scan.allxyz".format(base)) + angles, energies = self._read_orca_dat("{}_scan.relaxscanact.dat".format(base)) + if os.path.isfile("{}_sp.xyzact.dat".format(base)): + _, energies = self._read_orca_dat("{}_sp.xyzact.dat".format(base)) energies = np.array(energies) * Hartree * mol / kJ + return n_atoms, coords, angles, energies, point_charges