From 223ec89fc757614ad0ecd9892c78ca85dd53f786 Mon Sep 17 00:00:00 2001 From: ericwang6 Date: Tue, 27 Apr 2021 01:04:36 +0800 Subject: [PATCH 01/13] support read and dump multiframes for .gro file --- dpdata/gromacs/gro.py | 83 +++++++++++++++------ tests/gromacs/case_for_format_atom_name.gro | 35 +++++++++ tests/gromacs/multi_frames.gro | 24 ++++++ 3 files changed, 118 insertions(+), 24 deletions(-) create mode 100644 tests/gromacs/case_for_format_atom_name.gro create mode 100644 tests/gromacs/multi_frames.gro diff --git a/dpdata/gromacs/gro.py b/dpdata/gromacs/gro.py index 6e12fb1e1..35879f873 100644 --- a/dpdata/gromacs/gro.py +++ b/dpdata/gromacs/gro.py @@ -1,11 +1,23 @@ #!/usr/bin/env python3 - +import re import numpy as np +import json +import os nm2ang = 10. +ang2nm = 0.1 +cell_idx_gmx2dp = [0, 4, 8, 1, 2, 3, 5, 6, 7] + +def _format_atom_name(atom_name): + patt = re.compile("[a-zA-Z]*") + match = re.search(patt, atom_name) + fmt_name = match.group().capitalize() + return fmt_name -def _get_line(line): +def _get_line(line, fmt_atom_name=True): atom_name = line[10:15].split()[0] + if fmt_atom_name: + atom_name = _format_atom_name(atom_name) atom_idx = int(line[15:20].split()[0]) posis = [float(line[ii:ii+8]) for ii in range(20,44,8)] posis = np.array(posis) * nm2ang @@ -29,27 +41,50 @@ def _get_cell(line): cell = cell * nm2ang return cell -def file_to_system_data(fname): - names = [] - idxs = [] - posis = [] +def file_to_system_data(fname, format_atom_name=True): + system = {'coords': [], 'cells': []} with open(fname) as fp: - fp.readline() - natoms = int(fp.readline()) - for ii in range(natoms): - n, i, p = _get_line(fp.readline()) - names.append(n) - idxs.append(i) - posis.append(p) - cell = _get_cell(fp.readline()) - posis = np.array(posis) - system = {} - system['orig'] = np.array([0, 0, 0]) - system['atom_names'] = list(set(names)) - system['atom_names'].sort() - system['atom_numbs'] = [names.count(ii) for ii in system['atom_names']] - system['atom_types'] = [system['atom_names'].index(ii) for ii in names] - system['atom_types'] = np.array(system['atom_types'], dtype = int) - system['coords'] = np.array([posis]) - system['cells'] = np.array([cell]) + frame = 0 + while True: + flag = fp.readline() + if not flag: + break + else: + frame += 1 + names = [] + idxs = [] + posis = [] + natoms = int(fp.readline()) + for ii in range(natoms): + n, i, p = _get_line(fp.readline(), fmt_atom_name=format_atom_name) + names.append(n) + idxs.append(i) + posis.append(p) + cell = _get_cell(fp.readline()) + posis = np.array(posis) + if frame == 1: + system['orig'] = np.zeros(3) + system['atom_names'] = list(set(names)) + system['atom_numbs'] = [names.count(ii) for ii in system['atom_names']] + system['atom_types'] = [system['atom_names'].index(ii) for ii in names] + system['atom_types'] = np.array(system['atom_types'], dtype = int) + system['coords'].append(posis) + system['cells'].append(cell) + system['coords'] = np.array(system['coords']) + system['cells'] = np.array(system['cells']) return system + +def from_system_data(system, f_idx=0): + ret = "" + ret += " molecule" + "\n" + n_atoms = sum(system["atom_numbs"]) + ret += " " + str(n_atoms) + "\n" + for i in range(n_atoms): + atom_type = system["atom_types"][i] + atom_name = system["atom_names"][atom_type] + coords = system["coords"][f_idx] * ang2nm + ret += "{:>5d}{:<5s}{:>5s}{:5d}{:8.3f}{:8.3f}{:8.3f}\n".format(1, "MOL", atom_name, i+1, *tuple(coords[i])) + cell = (system["cells"][f_idx].flatten() * ang2nm)[cell_idx_gmx2dp] + ret += " " + " ".join([str(x) for x in cell]) + + return ret \ No newline at end of file diff --git a/tests/gromacs/case_for_format_atom_name.gro b/tests/gromacs/case_for_format_atom_name.gro new file mode 100644 index 000000000..29b248bd4 --- /dev/null +++ b/tests/gromacs/case_for_format_atom_name.gro @@ -0,0 +1,35 @@ +ligand in water t= 0.00000 step= 0 + 32 + 1MOL C1 1 2.763 2.029 0.966 + 1MOL C2 2 2.844 2.141 0.941 + 1MOL C3 3 2.788 2.268 0.958 + 1MOL C4 4 2.653 2.285 0.991 + 1MOL C5 5 2.576 2.171 1.011 + 1MOL C6 6 2.630 2.042 1.007 + 1MOL CL1 7 2.405 2.194 1.045 + 1MOL C7 8 2.588 2.419 0.999 + 1MOL O1 9 2.558 2.474 0.894 + 1MOL N1 10 2.572 2.471 1.125 + 1MOL C8 11 2.504 2.587 1.170 + 1MOL C9 12 2.542 2.648 1.289 + 1MOL C10 13 2.461 2.753 1.336 + 1MOL N2 14 2.359 2.804 1.266 + 1MOL C11 15 2.321 2.742 1.153 + 1MOL C12 16 2.389 2.631 1.103 + 1MOL N3 17 2.217 2.806 1.082 + 1MOL C13 18 2.118 2.742 1.011 + 1MOL O2 19 2.121 2.624 0.980 + 1MOL C14 20 2.003 2.827 0.960 + 1MOL CL2 21 2.875 2.415 0.920 + 1MOL H1 22 2.026 2.932 0.936 + 1MOL H2 23 2.810 1.932 0.951 + 1MOL H3 24 2.945 2.132 0.902 + 1MOL H4 25 2.577 1.954 1.041 + 1MOL H5 26 2.608 2.408 1.196 + 1MOL H6 27 2.624 2.605 1.346 + 1MOL H7 28 2.483 2.804 1.430 + 1MOL H8 29 2.357 2.582 1.012 + 1MOL H9 30 2.203 2.904 1.105 + 1MOL H10 31 1.920 2.816 1.031 + 1MOL H11 32 1.967 2.787 0.864 + 3.24290 3.24290 2.29308 0.00000 0.00000 0.00000 0.00000 1.62145 1.62145 \ No newline at end of file diff --git a/tests/gromacs/multi_frames.gro b/tests/gromacs/multi_frames.gro new file mode 100644 index 000000000..2f16e4d6a --- /dev/null +++ b/tests/gromacs/multi_frames.gro @@ -0,0 +1,24 @@ +frame_1 +9 + 1SOL O 1 0.135 0.183 0.341 + 1SOL H 2 0.177 0.149 0.262 + 1SOL H 3 0.046 0.149 0.339 + 2SOL O 4 0.520 0.447 0.111 + 2SOL H 5 0.567 0.481 0.035 + 2SOL H 6 0.568 0.481 0.186 + 3SOL O 7 0.651 0.539 0.335 + 3SOL H 8 0.653 0.634 0.336 + 3SOL H 9 0.743 0.512 0.336 + 0.7822838765564372 0.7353572647182051 0.9036518515423753 +frame_2 +9 + 1SOL O 1 0.135 0.183 0.341 + 1SOL H 2 0.177 0.149 0.262 + 1SOL H 3 0.046 0.149 0.339 + 2SOL O 4 0.520 0.447 0.111 + 2SOL H 5 0.567 0.481 0.035 + 2SOL H 6 0.568 0.481 0.186 + 3SOL O 7 0.651 0.539 0.335 + 3SOL H 8 0.653 0.634 0.336 + 3SOL H 9 0.743 0.512 0.336 + 0.7822838765564372 0.7353572647182051 0.9036518515423753 From c660136e9afc5070591e7cc6d9d847ef96609b80 Mon Sep 17 00:00:00 2001 From: ericwang6 Date: Tue, 27 Apr 2021 01:05:29 +0800 Subject: [PATCH 02/13] add unittest case for gromacs/gro.py --- tests/test_gromacs_gro.py | 75 ++++++++++++++++++++++++++++++++++++--- 1 file changed, 71 insertions(+), 4 deletions(-) diff --git a/tests/test_gromacs_gro.py b/tests/test_gromacs_gro.py index 1acb4d27f..16070448b 100644 --- a/tests/test_gromacs_gro.py +++ b/tests/test_gromacs_gro.py @@ -5,8 +5,9 @@ class TestGromacsGro(unittest.TestCase): def test_read_file(self): - system = dpdata.System('gromacs/1h.gro') - self.assertEqual(system['atom_names'], ['H', 'O']) + system = dpdata.System('gromacs/1h.gro', type_map=['H', 'O']) + self.assertTrue('H' in system['atom_names']) + self.assertTrue('O' in system['atom_names']) self.assertEqual(system['atom_numbs'], [6, 3]) for cc,ii in enumerate([1, 0, 0, 1, 0, 0, 1, 0, 0]): self.assertEqual(system['atom_types'][cc], ii) @@ -24,8 +25,9 @@ def test_read_file(self): self.assertAlmostEqual(system['coords'][0][8][2], 3.36) def test_read_file_tri(self): - system = dpdata.System('gromacs/1h.tri.gro') - self.assertEqual(system['atom_names'], ['H', 'O']) + system = dpdata.System('gromacs/1h.tri.gro', type_map=['H', 'O']) + self.assertTrue('H' in system['atom_names']) + self.assertTrue('O' in system['atom_names']) self.assertEqual(system['atom_numbs'], [6, 3]) for cc,ii in enumerate([1, 0, 0, 1, 0, 0, 1, 0, 0]): self.assertEqual(system['atom_types'][cc], ii) @@ -44,3 +46,68 @@ def test_read_file_tri(self): self.assertAlmostEqual(system['coords'][0][8][1], 5.12) self.assertAlmostEqual(system['coords'][0][8][2], 3.36) system.to('vasp/poscar', 'POSCAR') + +class TestGromacsGroMultiFrames(unittest.TestCase): + def test_read_file(self): + system = dpdata.System('gromacs/multi_frames.gro', type_map=['H', 'O']) + self.assertTrue('H' in system['atom_names']) + self.assertTrue('O' in system['atom_names']) + self.assertEqual(system['atom_numbs'], [6, 3]) + for cc,ii in enumerate([1, 0, 0, 1, 0, 0, 1, 0, 0]): + self.assertEqual(system['atom_types'][cc], ii) + self.assertEqual(len(system['cells']), 2) + self.assertEqual(len(system['coords']), 2) + for ii in range(3): + for jj in range(3): + if ii != jj: + self.assertAlmostEqual(system['cells'][0][ii][jj], 0) # frame no.1 + self.assertAlmostEqual(system['cells'][1][ii][jj], 0) # frame no.2 + # frame no.1 + self.assertAlmostEqual(system['cells'][0][0][0], 7.822838765564372) + self.assertAlmostEqual(system['cells'][0][1][1], 7.353572647182051) + self.assertAlmostEqual(system['cells'][0][2][2], 9.036518515423753) + self.assertAlmostEqual(system['coords'][0][8][0], 7.43) + self.assertAlmostEqual(system['coords'][0][8][1], 5.12) + self.assertAlmostEqual(system['coords'][0][8][2], 3.36) + # frame no.2 + self.assertAlmostEqual(system['cells'][1][0][0], 7.822838765564372) + self.assertAlmostEqual(system['cells'][1][1][1], 7.353572647182051) + self.assertAlmostEqual(system['cells'][1][2][2], 9.036518515423753) + self.assertAlmostEqual(system['coords'][1][8][0], 7.43) + self.assertAlmostEqual(system['coords'][1][8][1], 5.12) + self.assertAlmostEqual(system['coords'][1][8][2], 3.36) + + +class TestFormatAtomName(unittest.TestCase): + def test_format_atom_name(self): + system = dpdata.System("gromacs/case_for_format_atom_name.gro", fmt='gromacs/gro', type_map=['H','C','N','O','Cl']) + self.assertEqual(system.formula, "H11C14N3O2Cl2") + + def test_no_format_atom_name(self): + system = dpdata.System("gromacs/case_for_format_atom_name.gro", fmt='gromacs/gro', format_atom_name=False) + atoms = ['CL1', 'H6', 'C4', 'C3', 'C6', 'C11', 'H10', 'C2', 'N3', 'C14', + 'H7', 'H8', 'C13', 'H2', 'H1', 'H4', 'O2', 'H9', 'O1', 'N2', 'C9', + 'H3', 'C5', 'H11', 'N1', 'C7', 'C10', 'CL2', 'H5', 'C1', 'C8','C12'] + for at in atoms: + self.assertTrue(at in system['atom_names']) + + +class TestDumpGromacsGro(unittest.TestCase): + def setUp(self): + self.system = dpdata.System('gromacs/multi_frames.gro', type_map=['H', 'O']) + + def test_dump_single_frame(self): + self.system.to_gromacs_gro('gromacs/tmp_1.gro', frame_idx=0) + tmp = dpdata.System('gromacs/tmp_1.gro', type_map=['H', 'O']) + self.assertEqual(tmp.get_nframes(), 1) + + def test_dump_multi_frames(self): + self.system.to_gromacs_gro('gromacs/tmp_2.gro') + tmp = dpdata.System('gromacs/tmp_2.gro', type_map=['H', 'O']) + self.assertEqual(tmp.get_nframes(), 2) + + def tearDown(self): + if os.path.exists('gromacs/tmp_1.gro'): + os.remove('gromacs/tmp_1.gro') + if os.path.exists('gromacs/tmp_2.gro'): + os.remove('gromacs/tmp_2.gro') From 28313fa09a9813b589157acf70894cc44c98a3da Mon Sep 17 00:00:00 2001 From: ericwang6 Date: Tue, 27 Apr 2021 01:06:40 +0800 Subject: [PATCH 03/13] add dump system to .gjf file for gaussian --- dpdata/gaussian/gjf.py | 16 ++++++++ dpdata/system.py | 91 +++++++++++++++++++++++++++++++++++++++++- 2 files changed, 105 insertions(+), 2 deletions(-) create mode 100644 dpdata/gaussian/gjf.py diff --git a/dpdata/gaussian/gjf.py b/dpdata/gaussian/gjf.py new file mode 100644 index 000000000..1e8be587c --- /dev/null +++ b/dpdata/gaussian/gjf.py @@ -0,0 +1,16 @@ +def to_gjf_string(system, frame_idx=0, header="", title="", foot="", charge=0, mult=1): + if not title: + title = system.formula + if not header: + header = "#force B3LYP/6-31G(d)" + ret = header + "\n\n" + title + "\n\n" + str(charge) + " " + str(mult) + "\n" + coord = system.data["coords"][frame_idx].reshape(-1, 3) + for atype, c in zip(system.data["atom_types"], coord): + ret += system.data["atom_names"][atype] + " " + ret += " ".join([str(x) for x in c]) + ret += "\n" + if foot: + ret += "\n" + foot + "\n" + else: + ret += "\n" + return ret \ No newline at end of file diff --git a/dpdata/system.py b/dpdata/system.py index 0af601896..abb6c061d 100644 --- a/dpdata/system.py +++ b/dpdata/system.py @@ -15,6 +15,7 @@ import dpdata.siesta.output import dpdata.siesta.aiMD_output import dpdata.md.pbc +import dpdata.gaussian.gjf import dpdata.gaussian.log import dpdata.amber.md import dpdata.cp2k.output @@ -622,7 +623,7 @@ def from_deepmd_raw(self, folder, type_map = None) : @register_from_funcs.register_funcs("gro") @register_from_funcs.register_funcs("gromacs/gro") - def from_gromacs_gro(self, file_name) : + def from_gromacs_gro(self, file_name, format_atom_name=True) : """ Load gromacs .gro file @@ -631,7 +632,44 @@ def from_gromacs_gro(self, file_name) : file_name : str The input file name """ - self.data = dpdata.gromacs.gro.file_to_system_data(file_name) + self.data = dpdata.gromacs.gro.file_to_system_data(file_name, format_atom_name=format_atom_name) + + @register_to_funcs.register_funcs("gromacs/string") + def to_gromacs_string(self, frame_idx=-1): + """ + Dump the system in gromacs .gro format string + + Parameters + ---------- + frame_idx : int + The index of the frame to dump + """ + assert(frame_idx < len(self.data['coords'])) + if frame_idx == -1: + strs = [] + for idx in range(self.get_nframes()): + gro_str = dpdata.gromacs.gro.from_system_data(self.data, f_idx=idx) + strs.append(gro_str) + return "\n".join(strs) + else: + gro_str = dpdata.gromacs.gro.from_system_data(self.data, f_idx=frame_idx) + return gro_str + + @register_to_funcs.register_funcs("gromacs/gro") + def to_gromacs_gro(self, file_name, frame_idx=-1): + """ + Dump the system in gromacs .gro format + + Parameters + ---------- + file_name : str + The output file name + frame_idx : int + The index of the frame to dump + """ + gro_str=self.to_gromacs_string(frame_idx=frame_idx) + with open(file_name, 'w') as fp: + fp.write(gro_str) @register_to_funcs.register_funcs("deepmd/npy") def to_deepmd_npy(self, folder, set_size = 5000, prec=np.float32) : @@ -709,7 +747,56 @@ def to_pwmat_atomconfig(self, file_name, frame_idx = 0) : w_str = dpdata.pwmat.atomconfig.from_system_data(self.data, frame_idx) with open(file_name, 'w') as fp: fp.write(w_str) + + @register_to_funcs.register_funcs("gaussian/str") + def to_gaussian_str(self, frame_idx=0, header="", title="", foot="", charge=0, mult=1): + """ + Write input file for gaussian + + Parameters + ---------- + frame_idx: int + The index of the frame to specify molecule geometry + header: str + The gaussian input file headers, include link0 and route section + title: str + The title section of gaussian input file + foot: str + Additional section for 'Opt=ModRedundant' keyword + charge: int + The charge of molecule + mult: int + The multiplicity of molecule + """ + assert frame_idx < self.get_nframes() + ret = dpdata.gaussian.gjf.to_gjf_string(self, frame_idx, header, title, foot, charge, mult) + + return ret + + @register_to_funcs.register_funcs("gaussian/gjf") + def to_gaussian_gjf(self, gjf_file, frame_idx=0, header="", title="", foot="", charge=0, mult=1): + """ + Write input file for gaussian + Parameters + ---------- + gjf_file: str + The path of .gjf file to write + frame_idx: int + The index of the frame to specify molecule geometry + header: str + The gaussian input file headers, include link0 and route section + title: str + The title section of gaussian input file + foot: str + Additional section for 'Opt=ModRedundant' keyword + charge: int + The charge of molecule + mult: int + The multiplicity of molecule + """ + with open(gjf_file, 'w+') as f: + f.write(self.to_gaussian_str(frame_idx, header, title, foot, charge, mult)) def affine_map(self, trans, f_idx = 0) : assert(np.linalg.det(trans) != 0) From 0a2de6857ed01fd1cdf6a2ce6a6ac49e1cf02dea Mon Sep 17 00:00:00 2001 From: ericwang6 Date: Tue, 27 Apr 2021 01:07:25 +0800 Subject: [PATCH 04/13] add unittest case for write .gjf file --- tests/test_gaussian_gjf.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) create mode 100644 tests/test_gaussian_gjf.py diff --git a/tests/test_gaussian_gjf.py b/tests/test_gaussian_gjf.py new file mode 100644 index 000000000..b2c1737d3 --- /dev/null +++ b/tests/test_gaussian_gjf.py @@ -0,0 +1,16 @@ +import os +import unittest +from context import dpdata + + +class TestDumpGaussianGjf(unittest.TestCase): + def setUp(self): + self.system = dpdata.LabeledSystem('gaussian/methane.gaussianlog', + fmt = 'gaussian/log') + + def test_dump_to_gjf(self): + self.system.to_gaussian_gjf("gaussian/tmp.gjf") + + def tearDown(self): + if os.path.exists('gaussian/tmp.gjf'): + os.remove('gaussian/tmp.gjf') \ No newline at end of file From cca873d35c6e2bed3012a047ae2b4ae89f54fb16 Mon Sep 17 00:00:00 2001 From: ericwang6 Date: Tue, 27 Apr 2021 01:07:51 +0800 Subject: [PATCH 05/13] support multi-frames for .gro file --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 8dd9c8649..1f4b33387 100644 --- a/README.md +++ b/README.md @@ -78,7 +78,7 @@ The `System` or `LabeledSystem` can be constructed from the following file forma | PWmat | movement | True | True | LabeledSystem | 'pwmat/movement' | | PWmat | OUT.MLMD | True | True | LabeledSystem | 'pwmat/out.mlmd' | | Amber | multi | True | True | LabeledSystem | 'amber/md' | -| Gromacs | gro | False | False | System | 'gromacs/gro' | +| Gromacs | gro | True | False | System | 'gromacs/gro' | The Class `dpdata.MultiSystems` can read data from a dir which may contains many files of different systems, or from single xyz file which contains different systems. From 2435d63386adcf2000b513d25f507122ac533900 Mon Sep 17 00:00:00 2001 From: Ericwang6 <63655850+Ericwang6@users.noreply.github.com> Date: Tue, 27 Apr 2021 10:50:20 +0800 Subject: [PATCH 06/13] change ang2nm = 0.1 to ang2nm = 1 / nm2ang --- dpdata/gromacs/gro.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dpdata/gromacs/gro.py b/dpdata/gromacs/gro.py index 35879f873..90cc4096d 100644 --- a/dpdata/gromacs/gro.py +++ b/dpdata/gromacs/gro.py @@ -5,7 +5,7 @@ import os nm2ang = 10. -ang2nm = 0.1 +ang2nm = 1 / nm2ang cell_idx_gmx2dp = [0, 4, 8, 1, 2, 3, 5, 6, 7] def _format_atom_name(atom_name): @@ -87,4 +87,4 @@ def from_system_data(system, f_idx=0): cell = (system["cells"][f_idx].flatten() * ang2nm)[cell_idx_gmx2dp] ret += " " + " ".join([str(x) for x in cell]) - return ret \ No newline at end of file + return ret From ff46250b0dc6f2cc466696073ed2ef883141965e Mon Sep 17 00:00:00 2001 From: Ericwang6 <63655850+Ericwang6@users.noreply.github.com> Date: Tue, 27 Apr 2021 10:51:59 +0800 Subject: [PATCH 07/13] Update gro.py --- dpdata/gromacs/gro.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dpdata/gromacs/gro.py b/dpdata/gromacs/gro.py index 90cc4096d..b2473a03b 100644 --- a/dpdata/gromacs/gro.py +++ b/dpdata/gromacs/gro.py @@ -5,7 +5,7 @@ import os nm2ang = 10. -ang2nm = 1 / nm2ang +ang2nm = 1. / nm2ang cell_idx_gmx2dp = [0, 4, 8, 1, 2, 3, 5, 6, 7] def _format_atom_name(atom_name): From 32993e48d1c6e7232a19a1f046145c6325673a4e Mon Sep 17 00:00:00 2001 From: Ericwang6 Date: Tue, 27 Apr 2021 11:26:05 +0800 Subject: [PATCH 08/13] modify ang2nm --- dpdata/gromacs/gro.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dpdata/gromacs/gro.py b/dpdata/gromacs/gro.py index 35879f873..081851d3f 100644 --- a/dpdata/gromacs/gro.py +++ b/dpdata/gromacs/gro.py @@ -5,7 +5,7 @@ import os nm2ang = 10. -ang2nm = 0.1 +ang2nm = 1. / nm2ang cell_idx_gmx2dp = [0, 4, 8, 1, 2, 3, 5, 6, 7] def _format_atom_name(atom_name): From c4796fcdf6eb5ff4ae1ac055da544e63c4a44825 Mon Sep 17 00:00:00 2001 From: Ericwang6 Date: Tue, 27 Apr 2021 11:27:23 +0800 Subject: [PATCH 09/13] add doc string --- dpdata/gaussian/gjf.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/dpdata/gaussian/gjf.py b/dpdata/gaussian/gjf.py index 1e8be587c..0c3641034 100644 --- a/dpdata/gaussian/gjf.py +++ b/dpdata/gaussian/gjf.py @@ -1,4 +1,24 @@ def to_gjf_string(system, frame_idx=0, header="", title="", foot="", charge=0, mult=1): + ''' + Generate a string that can be written to a Gaussian input file, please refer to https://gaussian.com/input/?tabid=2 for more information on syntax of Gaussian input files. + + Parameters + -------- + system : dpdata.System + The system to write a gaussian input file + frame_idx : int + The index of frame to specify molecule geometry + header : str + The route section (# lines) and link-0 commands (% commands) + title : str + The title section of .gjf files + foot : str + The modifications to coordinates, used when Opt=ModRedundant + charge : int + The charge of the system + mult : int + The multiplicity of the system + ''' if not title: title = system.formula if not header: From b75bb3d4e39c012747799168dfc62f5dac50b7f1 Mon Sep 17 00:00:00 2001 From: Ericwang6 Date: Tue, 27 Apr 2021 11:28:29 +0800 Subject: [PATCH 10/13] merge gromacs/string, gaussian/string to gromacs/gro, gaussian/gjf --- dpdata/system.py | 96 +++++++++++++++++------------------------------- 1 file changed, 34 insertions(+), 62 deletions(-) diff --git a/dpdata/system.py b/dpdata/system.py index abb6c061d..763bc5bd8 100644 --- a/dpdata/system.py +++ b/dpdata/system.py @@ -634,13 +634,15 @@ def from_gromacs_gro(self, file_name, format_atom_name=True) : """ self.data = dpdata.gromacs.gro.file_to_system_data(file_name, format_atom_name=format_atom_name) - @register_to_funcs.register_funcs("gromacs/string") - def to_gromacs_string(self, frame_idx=-1): + @register_to_funcs.register_funcs("gromacs/gro") + def to_gromacs_gro(self, file_name=None, frame_idx=-1): """ - Dump the system in gromacs .gro format string - + Dump the system in gromacs .gro format + Parameters ---------- + file_name : str or None + The output file name. If None, return the file content as a string frame_idx : int The index of the frame to dump """ @@ -650,26 +652,15 @@ def to_gromacs_string(self, frame_idx=-1): for idx in range(self.get_nframes()): gro_str = dpdata.gromacs.gro.from_system_data(self.data, f_idx=idx) strs.append(gro_str) - return "\n".join(strs) + gro_str = "\n".join(strs) else: gro_str = dpdata.gromacs.gro.from_system_data(self.data, f_idx=frame_idx) + + if file_name is None: return gro_str - - @register_to_funcs.register_funcs("gromacs/gro") - def to_gromacs_gro(self, file_name, frame_idx=-1): - """ - Dump the system in gromacs .gro format - - Parameters - ---------- - file_name : str - The output file name - frame_idx : int - The index of the frame to dump - """ - gro_str=self.to_gromacs_string(frame_idx=frame_idx) - with open(file_name, 'w') as fp: - fp.write(gro_str) + else: + with open(file_name, 'w+') as fp: + fp.write(gro_str) @register_to_funcs.register_funcs("deepmd/npy") def to_deepmd_npy(self, folder, set_size = 5000, prec=np.float32) : @@ -748,55 +739,36 @@ def to_pwmat_atomconfig(self, file_name, frame_idx = 0) : with open(file_name, 'w') as fp: fp.write(w_str) - @register_to_funcs.register_funcs("gaussian/str") - def to_gaussian_str(self, frame_idx=0, header="", title="", foot="", charge=0, mult=1): + @register_to_funcs.register_funcs("gaussian/gjf") + def to_gaussian_gjf(self, gjf_file=None, frame_idx=0, header="", title="", foot="", charge=0, mult=1): """ - Write input file for gaussian + Write input file for Gaussian. Please refer to https://gaussian.com/input/?tabid=2 for more information on syntax of Gaussian input files. Parameters ---------- - frame_idx: int - The index of the frame to specify molecule geometry - header: str - The gaussian input file headers, include link0 and route section - title: str - The title section of gaussian input file - foot: str - Additional section for 'Opt=ModRedundant' keyword - charge: int - The charge of molecule - mult: int - The multiplicity of molecule + gjf_file : str or None + The file to write. If None, return the file's content as a string + frame_idx : int + The index of frame to specify molecule geometry + header : str + The route section (# lines) and link-0 commands (% commands) + title : str + The title section of .gjf files + foot : str + The modifications to coordinates, used when Opt=ModRedundant + charge : int + The charge of the system + mult : int + The multiplicity of the system """ assert frame_idx < self.get_nframes() ret = dpdata.gaussian.gjf.to_gjf_string(self, frame_idx, header, title, foot, charge, mult) - return ret - - @register_to_funcs.register_funcs("gaussian/gjf") - def to_gaussian_gjf(self, gjf_file, frame_idx=0, header="", title="", foot="", charge=0, mult=1): - """ - Write input file for gaussian - - Parameters - ---------- - gjf_file: str - The path of .gjf file to write - frame_idx: int - The index of the frame to specify molecule geometry - header: str - The gaussian input file headers, include link0 and route section - title: str - The title section of gaussian input file - foot: str - Additional section for 'Opt=ModRedundant' keyword - charge: int - The charge of molecule - mult: int - The multiplicity of molecule - """ - with open(gjf_file, 'w+') as f: - f.write(self.to_gaussian_str(frame_idx, header, title, foot, charge, mult)) + if gjf_file is None: + return ret + else: + with open(gjf_file, 'w+') as f: + f.write(ret) def affine_map(self, trans, f_idx = 0) : assert(np.linalg.det(trans) != 0) From 2d2d43db4b5069734aef78095424212a5105274e Mon Sep 17 00:00:00 2001 From: Ericwang6 Date: Wed, 28 Apr 2021 00:00:28 +0800 Subject: [PATCH 11/13] add test for write gjf files --- tests/test_gaussian_gjf.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/tests/test_gaussian_gjf.py b/tests/test_gaussian_gjf.py index b2c1737d3..08a0fae00 100644 --- a/tests/test_gaussian_gjf.py +++ b/tests/test_gaussian_gjf.py @@ -10,6 +10,27 @@ def setUp(self): def test_dump_to_gjf(self): self.system.to_gaussian_gjf("gaussian/tmp.gjf") + with open("gaussian/tmp.gjf") as f: + header = f.readline().strip() + f.readline() + title = f.readline().strip() + f.readline() + charge, mult = (int(x) for x in f.readline().strip().split()) + atoms = [] + coords = [] + for ii in range(5): + line = f.readline().strip().split() + atoms.append(line[0]) + coords.append([float(x) for x in line[1:]]) + + self.assertEqual(header, "#force B3LYP/6-31G(d)") + self.assertEqual(title, self.system.formula) + self.assertEqual(charge, 0) + self.assertEqual(mult, 1) + self.assertEqual(atoms, ['C', 'H', 'H', 'H', 'H']) + for i in range(self.system['coords'].shape[1]): + for j in range(3): + self.assertAlmostEqual(coords[i][j], self.system['coords'][0][i][j]) def tearDown(self): if os.path.exists('gaussian/tmp.gjf'): From c806a0a04383858a6d94836d023d626af56aa392 Mon Sep 17 00:00:00 2001 From: Ericwang6 Date: Thu, 29 Apr 2021 20:40:11 +0800 Subject: [PATCH 12/13] support read and dump multiframes for .gro file --- dpdata/system.py | 31 ------------------------------- 1 file changed, 31 deletions(-) diff --git a/dpdata/system.py b/dpdata/system.py index 763bc5bd8..5db810396 100644 --- a/dpdata/system.py +++ b/dpdata/system.py @@ -15,7 +15,6 @@ import dpdata.siesta.output import dpdata.siesta.aiMD_output import dpdata.md.pbc -import dpdata.gaussian.gjf import dpdata.gaussian.log import dpdata.amber.md import dpdata.cp2k.output @@ -738,37 +737,7 @@ def to_pwmat_atomconfig(self, file_name, frame_idx = 0) : w_str = dpdata.pwmat.atomconfig.from_system_data(self.data, frame_idx) with open(file_name, 'w') as fp: fp.write(w_str) - - @register_to_funcs.register_funcs("gaussian/gjf") - def to_gaussian_gjf(self, gjf_file=None, frame_idx=0, header="", title="", foot="", charge=0, mult=1): - """ - Write input file for Gaussian. Please refer to https://gaussian.com/input/?tabid=2 for more information on syntax of Gaussian input files. - Parameters - ---------- - gjf_file : str or None - The file to write. If None, return the file's content as a string - frame_idx : int - The index of frame to specify molecule geometry - header : str - The route section (# lines) and link-0 commands (% commands) - title : str - The title section of .gjf files - foot : str - The modifications to coordinates, used when Opt=ModRedundant - charge : int - The charge of the system - mult : int - The multiplicity of the system - """ - assert frame_idx < self.get_nframes() - ret = dpdata.gaussian.gjf.to_gjf_string(self, frame_idx, header, title, foot, charge, mult) - - if gjf_file is None: - return ret - else: - with open(gjf_file, 'w+') as f: - f.write(ret) def affine_map(self, trans, f_idx = 0) : assert(np.linalg.det(trans) != 0) From 471de612efdb38f2e09a16d7ee0e9956b6da1791 Mon Sep 17 00:00:00 2001 From: Ericwang6 Date: Thu, 29 Apr 2021 20:42:55 +0800 Subject: [PATCH 13/13] remove not fully developed support for .gjf file --- dpdata/gaussian/gjf.py | 36 ------------------------------------ tests/test_gaussian_gjf.py | 37 ------------------------------------- 2 files changed, 73 deletions(-) delete mode 100644 dpdata/gaussian/gjf.py delete mode 100644 tests/test_gaussian_gjf.py diff --git a/dpdata/gaussian/gjf.py b/dpdata/gaussian/gjf.py deleted file mode 100644 index 0c3641034..000000000 --- a/dpdata/gaussian/gjf.py +++ /dev/null @@ -1,36 +0,0 @@ -def to_gjf_string(system, frame_idx=0, header="", title="", foot="", charge=0, mult=1): - ''' - Generate a string that can be written to a Gaussian input file, please refer to https://gaussian.com/input/?tabid=2 for more information on syntax of Gaussian input files. - - Parameters - -------- - system : dpdata.System - The system to write a gaussian input file - frame_idx : int - The index of frame to specify molecule geometry - header : str - The route section (# lines) and link-0 commands (% commands) - title : str - The title section of .gjf files - foot : str - The modifications to coordinates, used when Opt=ModRedundant - charge : int - The charge of the system - mult : int - The multiplicity of the system - ''' - if not title: - title = system.formula - if not header: - header = "#force B3LYP/6-31G(d)" - ret = header + "\n\n" + title + "\n\n" + str(charge) + " " + str(mult) + "\n" - coord = system.data["coords"][frame_idx].reshape(-1, 3) - for atype, c in zip(system.data["atom_types"], coord): - ret += system.data["atom_names"][atype] + " " - ret += " ".join([str(x) for x in c]) - ret += "\n" - if foot: - ret += "\n" + foot + "\n" - else: - ret += "\n" - return ret \ No newline at end of file diff --git a/tests/test_gaussian_gjf.py b/tests/test_gaussian_gjf.py deleted file mode 100644 index 08a0fae00..000000000 --- a/tests/test_gaussian_gjf.py +++ /dev/null @@ -1,37 +0,0 @@ -import os -import unittest -from context import dpdata - - -class TestDumpGaussianGjf(unittest.TestCase): - def setUp(self): - self.system = dpdata.LabeledSystem('gaussian/methane.gaussianlog', - fmt = 'gaussian/log') - - def test_dump_to_gjf(self): - self.system.to_gaussian_gjf("gaussian/tmp.gjf") - with open("gaussian/tmp.gjf") as f: - header = f.readline().strip() - f.readline() - title = f.readline().strip() - f.readline() - charge, mult = (int(x) for x in f.readline().strip().split()) - atoms = [] - coords = [] - for ii in range(5): - line = f.readline().strip().split() - atoms.append(line[0]) - coords.append([float(x) for x in line[1:]]) - - self.assertEqual(header, "#force B3LYP/6-31G(d)") - self.assertEqual(title, self.system.formula) - self.assertEqual(charge, 0) - self.assertEqual(mult, 1) - self.assertEqual(atoms, ['C', 'H', 'H', 'H', 'H']) - for i in range(self.system['coords'].shape[1]): - for j in range(3): - self.assertAlmostEqual(coords[i][j], self.system['coords'][0][i][j]) - - def tearDown(self): - if os.path.exists('gaussian/tmp.gjf'): - os.remove('gaussian/tmp.gjf') \ No newline at end of file