diff --git a/dpdata/gaussian/gjf.py b/dpdata/gaussian/gjf.py index 727732b86..c2652747e 100644 --- a/dpdata/gaussian/gjf.py +++ b/dpdata/gaussian/gjf.py @@ -229,6 +229,12 @@ def make_gaussian_input( (symbol, frag_index[ii] + 1, *coordinate)) else: buff.append("%s %f %f %f" % (symbol, *coordinate)) + if not sys_data.get('nopbc', False): + # PBC condition + cell = sys_data['cells'][0] + for ii in range(3): + # use TV as atomic symbol, see https://gaussian.com/pbc/ + buff.append('TV %f %f %f' % (symbol, *cell[ii])) if basis_set is not None: # custom basis set buff.extend(['', basis_set, '']) diff --git a/dpdata/gaussian/log.py b/dpdata/gaussian/log.py index 5ed976fea..0e0066821 100644 --- a/dpdata/gaussian/log.py +++ b/dpdata/gaussian/log.py @@ -16,6 +16,8 @@ def to_system_data(file_name, md=False): coords_t = [] atom_symbols = [] forces_t = [] + cells_t = [] + nopbc = True with open(file_name) as fp: for line in fp: @@ -29,6 +31,7 @@ def to_system_data(file_name, md=False): flag = 5 coords = [] atom_symbols = [] + cells = [] if 1 <= flag <= 3 or 5 <= flag <= 9: flag += 1 @@ -38,18 +41,31 @@ def to_system_data(file_name, md=False): forces_t.append(forces) energy_t.append(energy) coords_t.append(coords) + if cells: + nopbc = False + cells_t.append(cells) + else: + cells_t.append([[100., 0., 0.], [0., 100., 0.], [0., 0., 100.]]) flag = 0 else: s = line.split() - forces.append([float(line[23:38]), float(line[38:53]), float(line[53:68])]) + if line[14:16] == "-2": + # PBC + pass + else: + forces.append([float(line[23:38]), float(line[38:53]), float(line[53:68])]) elif flag == 10: # atom_symbols and coords if line.startswith(" -------"): flag = 0 else: s = line.split() - coords.append([float(x) for x in s[3:6]]) - atom_symbols.append(symbols[int(s[1])]) + if int(s[1]) == -2: + # PBC cells, see https://gaussian.com/pbc/ + cells.append([float(x) for x in s[3:6]]) + else: + coords.append([float(x) for x in s[3:6]]) + atom_symbols.append(symbols[int(s[1])]) assert(coords_t), "cannot find coords" assert(energy_t), "cannot find energies" @@ -62,10 +78,11 @@ def to_system_data(file_name, md=False): forces_t = forces_t[-1:] energy_t = energy_t[-1:] coords_t = coords_t[-1:] + cells_t = cells_t[-1:] data['forces'] = np.array(forces_t) * force_convert data['energies'] = np.array(energy_t) * energy_convert data['coords'] = np.array(coords_t) data['orig'] = np.array([0, 0, 0]) - data['cells'] = np.array([[[100., 0., 0.], [0., 100., 0.], [0., 0., 100.]] for _ in energy_t]) - data['nopbc'] = True + data['cells'] = np.array(cells_t) + data['nopbc'] = nopbc return data diff --git a/tests/gaussian/h2pbc.gaussianlog b/tests/gaussian/h2pbc.gaussianlog new file mode 100644 index 000000000..b876e3171 --- /dev/null +++ b/tests/gaussian/h2pbc.gaussianlog @@ -0,0 +1,301 @@ + Entering Gaussian System, Link 0=g16 + Input=h.gjf + Output=h.log + Initial command: + /home/jz748/g16/g16/l1.exe "/home/jz748/tmp/Gau-1950599.inp" -scrdir="/home/jz748/tmp/" + Entering Link 1 = /home/jz748/g16/g16/l1.exe PID= 1950600. + + Copyright (c) 1988,1990,1992,1993,1995,1998,2003,2009,2016, + Gaussian, Inc. All Rights Reserved. + + This is part of the Gaussian(R) 16 program. It is based on + the Gaussian(R) 09 system (copyright 2009, Gaussian, Inc.), + the Gaussian(R) 03 system (copyright 2003, Gaussian, Inc.), + the Gaussian(R) 98 system (copyright 1998, Gaussian, Inc.), + the Gaussian(R) 94 system (copyright 1995, Gaussian, Inc.), + the Gaussian 92(TM) system (copyright 1992, Gaussian, Inc.), + the Gaussian 90(TM) system (copyright 1990, Gaussian, Inc.), + the Gaussian 88(TM) system (copyright 1988, Gaussian, Inc.), + the Gaussian 86(TM) system (copyright 1986, Carnegie Mellon + University), and the Gaussian 82(TM) system (copyright 1983, + Carnegie Mellon University). Gaussian is a federally registered + trademark of Gaussian, Inc. + + This software contains proprietary and confidential information, + including trade secrets, belonging to Gaussian, Inc. + + This software is provided under written license and may be + used, copied, transmitted, or stored only in accord with that + written license. + + The following legend is applicable only to US Government + contracts under FAR: + + RESTRICTED RIGHTS LEGEND + + Use, reproduction and disclosure by the US Government is + subject to restrictions as set forth in subparagraphs (a) + and (c) of the Commercial Computer Software - Restricted + Rights clause in FAR 52.227-19. + + Gaussian, Inc. + 340 Quinnipiac St., Bldg. 40, Wallingford CT 06492 + + + --------------------------------------------------------------- + Warning -- This program may not be used in any manner that + competes with the business of Gaussian, Inc. or will provide + assistance to any competitor of Gaussian, Inc. The licensee + of this program is prohibited from giving any competitor of + Gaussian, Inc. access to this program. By using this program, + the user acknowledges that Gaussian, Inc. is engaged in the + business of creating and licensing software in the field of + computational chemistry and represents and warrants to the + licensee that it is not a competitor of Gaussian, Inc. and that + it will not use this program in any manner prohibited above. + --------------------------------------------------------------- + + + Cite this work as: + Gaussian 16, Revision A.03, + M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, + M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, + G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, + J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, + J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, + F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, + T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, + G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, + J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, + T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, + F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, + V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, + K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, + J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, + J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, + J. B. Foresman, and D. J. Fox, Gaussian, Inc., Wallingford CT, 2016. + + ****************************************** + Gaussian 16: ES64L-G16RevA.03 25-Dec-2016 + 19-Dec-2022 + ****************************************** + -------------------- + # force wb97x/6-31g* + -------------------- + 1/10=7,30=1,38=1/1,3; + 2/12=2,17=6,18=5,40=1/2; + 3/5=1,6=6,7=1,11=2,25=1,30=1,71=1,74=-57/1,2,3; + 4//1; + 5/5=2,38=5/2; + 6/7=2,8=2,9=2,10=2,28=1/1; + 7/29=1/1,2,3,16; + 1/10=7,30=1/3; + 99//99; + - + H + - + Symbolic Z-matrix: + Charge = 0 Multiplicity = 1 + H 0. 0. 0. + H 1. 1. 0. + TV 10. 0. 0. + TV 0. 10. 0. + TV 0. 0. 10. + + + GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad + Berny optimization. + Initialization pass. + Trust Radius=3.00D-01 FncErr=1.00D-07 GrdErr=1.00D-07 EigMax=2.50D+02 EigMin=1.00D-04 + Number of steps in this run= 2 maximum allowed number of steps= 2. + GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad + + Before rotation: + --------------------------------------------------------------------- + Center Atomic Atomic Coordinates (Angstroms) + Number Number Type X Y Z + --------------------------------------------------------------------- + 1 1 0 0.000000 0.000000 0.000000 + 2 1 0 1.000000 1.000000 0.000000 + 3 -2 0 10.000000 0.000000 0.000000 + 4 -2 0 0.000000 10.000000 0.000000 + 5 -2 0 0.000000 0.000000 10.000000 + --------------------------------------------------------------------- + Lengths of translation vectors: 10.000000 10.000000 10.000000 + Angles of translation vectors: 90.000000 90.000000 90.000000 + --------------------------------------------------------------------- + Input orientation: + --------------------------------------------------------------------- + Center Atomic Atomic Coordinates (Angstroms) + Number Number Type X Y Z + --------------------------------------------------------------------- + 1 1 0 -0.500000 -0.500000 0.000000 + 2 1 0 0.500000 0.500000 0.000000 + 3 -2 0 10.000000 0.000000 0.000000 + 4 -2 0 0.000000 10.000000 0.000000 + 5 -2 0 0.000000 0.000000 10.000000 + --------------------------------------------------------------------- + Lengths of translation vectors: 10.000000 10.000000 10.000000 + Angles of translation vectors: 90.000000 90.000000 90.000000 + --------------------------------------------------------------------- + Distance matrix (angstroms): + 1 2 3 4 5 + 1 H 0.000000 + 2 H 1.414214 0.000000 + 3 TV 10.511898 9.513149 0.000000 + 4 TV 10.511898 9.513149 14.142136 0.000000 + 5 TV 10.024969 10.024969 14.142136 14.142136 0.000000 + Unit Cell Distance matrix (angstroms): + 1 2 + 1 H 0.000000 + 2 H 12.727922 0.000000 + Symmetry turned off: + Cannot cope with ghost atoms or with translation vectors. + Stoichiometry H2 + Framework group C1[X(H2)] + Deg. of freedom 0 + Full point group C1 NOp 1 + Standard basis: 6-31G(d) (6D, 7F) + 4 basis functions, 8 primitive gaussians, 4 cartesian basis functions + 1 alpha electrons 1 beta electrons + nuclear repulsion energy 0.3741847943 Hartrees. + NAtoms= 2 NActive= 2 NUniq= 2 SFac= 1.00D+00 NAtFMM= 60 NAOKFM=F Big=F + Integral buffers will be 131072 words long. + Raffenetti 2 integral format. + Two-electron integral symmetry is turned off. + FOutLm= 100.00. + Periodicity: 1 1 1 + Max integer dimensions: 6 6 6 + PBC vector 1 X= 18.8973 Y= 0.0000 Z= 0.0000 + PBC vector 2 X= 0.0000 Y= 18.8973 Z= 0.0000 + PBC vector 3 X= 0.0000 Y= 0.0000 Z= 18.8973 + Recp vector 1 X= 0.0529 Y= 0.0000 Z= 0.0000 + Recp vector 2 X= 0.0000 Y= 0.0529 Z= 0.0000 + Recp vector 3 X= 0.0000 Y= 0.0000 Z= 0.0529 + Generated k point mesh (from -Pi to Pi): + K space mesh: X= 14 Y= 14 Z= 14 + A half-cell shift: 0 + Using k point mesh (from -Pi to Pi): + K space mesh: X= 14 Y= 14 Z= 14 + A half-cell shift: 0 + CountK=T Total number of k points: 0 + CountK=T Total number of k points: 1376 + One-electron integrals computed using PRISM. + NBasis= 4 RedAO= T EigKep= 2.13D-01 NBF= 4 + NBsUse= 4 1.00D-06 EigRej= -1.00D+00 NBFU= 4 + RepCel: MaxNCR= 37 NClRep= 37 NMtPBC= 1389. + Harris functional with IExCor= 4538 and IRadAn= 5 diagonalized for initial guess. + RepCel: MaxNCR= 37 NClRep= 37 NMtPBC= 1389. + HarFok: IExCor= 4538 AccDes= 0.00D+00 IRadAn= 5 IDoV= 2 UseB2=F ITyADJ=14 + ICtDFT= 3500011 ScaDFX= 1.000000 1.000000 1.000000 1.000000 + FoFCou: FMM=T IPFlag= 524288 FMFlag= 990000 FMFlg1= 1001 + NFxFlg= 0 DoJE=T BraDBF=F KetDBF=T FulRan=T + wScrn= 0.000000 ICntrl= 500 IOpCl= 0 I1Cent= 200000004 NGrid= 0 + NMat0= 1 NMatS0= 1 NMatT0= 0 NMatD0= 1 NMtDS0= 0 NMtDT0= 0 + Symmetry not used in FoFCou. + FMM levels: 3 Number of levels for PrismC: 2 + Requested convergence on RMS density matrix=1.00D-07 within 128 cycles. + Requested convergence on MAX density matrix=1.00D-05. + Requested convergence on energy=1.00D-05. + No special actions if energy rises. + Diagonalized old Fock matrix for initial guess. + RepCel: MaxNCR= 37 NClRep= 37 NMtPBC= 1389. + RepCel: MaxNCR= 37 NClRep= 37 NMtPBC= 1389. + RepCel: MaxNCR= 37 NClRep= 37 NMtPBC= 1389. + RepCel: MaxNCR= 37 NClRep= 37 NMtPBC= 1389. + RepCel: MaxNCR= 37 NClRep= 37 NMtPBC= 1389. + RepCel: MaxNCR= 37 NClRep= 37 NMtPBC= 1389. + SCF Done: E(RwB97X) = -1.06463933888 A.U. after 6 cycles + NFock= 6 Conv=0.16D-08 -V/T= 2.3027 + + ********************************************************************** + + Population analysis using the SCF density. + + ********************************************************************** + + Condensed to atoms (all electrons): + 1 2 + 1 H 0.713229 0.286771 + 2 H 0.286771 0.713229 + Mulliken charges: + 1 + 1 H 0.000000 + 2 H -0.000000 + Sum of Mulliken charges = -0.00000 + Mulliken charges with hydrogens summed into heavy atoms: + 1 + Calling FoFJK, ICntrl= 2127 FMM=T ISym2X=0 I1Cent= 0 IOpClX= 0 NMat=1 NMatS=1 NMatT=0. + RepCel: MaxNCR= 37 NClRep= 37 NMtPBC= 1389. + ------------------------------------------------------------------- + Center Atomic Forces (Hartrees/Bohr) + Number Number X Y Z + ------------------------------------------------------------------- + 1 1 0.066019420 0.066019420 0.000000000 + 2 1 -0.066019420 -0.066019420 -0.000000000 + -2 -0.000000162 -0.000000016 0.000000000 + -2 -0.000000016 -0.000000162 -0.000000000 + -2 0.000000000 0.000000000 0.000000123 + ------------------------------------------------------------------- + Cartesian Forces: Max 0.066019420 RMS 0.034092282 + ----------------------------------------------------------------------------------------------- + Internal Coordinate Forces (Hartree/Bohr or radian) + Cent Atom N1 Length/X N2 Alpha/Y N3 Beta/Z J + ----------------------------------------------------------------------------------------------- + 1 H 0.066019( 1) 0.066019( 6) 0.000000( 11) + 2 H -0.066019( 2) -0.066019( 7) -0.000000( 12) + TV -0.000000( 3) -0.000000( 8) 0.000000( 13) + TV -0.000000( 4) -0.000000( 9) -0.000000( 14) + TV 0.000000( 5) 0.000000( 10) 0.000000( 15) + ----------------------------------------------------------------------------------------------- + Internal Forces: Max 0.066019420 RMS 0.034092282 + + GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad + Berny optimization. + Search for a local minimum. + Step number 1 out of a maximum of 2 + All quantities printed in internal units (Hartrees-Bohrs-Radians) + Second derivative matrix not updated -- first step. + ITU= 0 + Angle between quadratic step and forces= 0.00 degrees. + Linear search not attempted -- first point. + Variable Old X -DE/DX Delta X Delta X Delta X New X + (Linear) (Quad) (Total) + X1 -0.94486 0.06602 0.00000 0.06602 0.06602 -0.87884 + Y1 -0.94486 0.06602 0.00000 0.06602 0.06602 -0.87884 + Z1 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + X2 0.94486 -0.06602 0.00000 -0.06602 -0.06602 0.87884 + Y2 0.94486 -0.06602 0.00000 -0.06602 -0.06602 0.87884 + Z2 0.00000 -0.00000 0.00000 -0.00000 -0.00000 -0.00000 + X3 18.89726 -0.00000 0.00000 -0.00000 -0.00000 18.89726 + Y3 0.00000 -0.00000 0.00000 -0.00000 -0.00000 -0.00000 + Z3 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + X4 0.00000 -0.00000 0.00000 -0.00000 -0.00000 -0.00000 + Y4 18.89726 -0.00000 0.00000 -0.00000 -0.00000 18.89726 + Z4 0.00000 -0.00000 0.00000 -0.00000 -0.00000 -0.00000 + X5 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + Y5 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + Z5 18.89726 0.00000 0.00000 0.00000 0.00000 18.89726 + Item Value Threshold Converged? + Maximum Force 0.066019 0.000450 NO + RMS Force 0.034092 0.000300 NO + Maximum Displacement 0.066019 0.001800 NO + RMS Displacement 0.034092 0.001200 NO + Predicted change in Energy=-8.717128D-03 + GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad + + 1\1\GINC-LOCALHOST\Force\RwB97X\6-31G(d)\H2\JZ748\19-Dec-2022\0\\# for + ce wb97x/6-31g*\\H\\0,1\H,-0.5,-0.5,0.\H,0.5,0.5,0.\TV,10.,0.,0.\TV,0. + ,10.,0.\TV,0.,0.,10.\\Version=ES64L-G16RevA.03\HF=-1.0646393\RMSD=1.63 + 1e-09\RMSF=3.409e-02\PG=C01 [X(H2)]\\@ + + + ABOVE ALL I AM AN OPTIMIST FOR NUMBER THEORY, + AND I HOLD THE HOPE THAT WE MAY NOT BE FAR FROM + A TIME WHEN IRREFUTABLE ARITHMETIC WILL CELEBRATE + ITS TRIUMPHS IN PHYSICS AND CHEMISTRY. + -- HERMANN MINKOWSKI, 1905 + Job cpu time: 0 days 0 hours 0 minutes 2.4 seconds. + Elapsed time: 0 days 0 hours 0 minutes 2.3 seconds. + File lengths (MBytes): RWF= 16 Int= 0 D2E= 0 Chk= 2 Scr= 1 + Normal termination of Gaussian 16 at Mon Dec 19 01:52:07 2022. diff --git a/tests/test_gaussian_log.py b/tests/test_gaussian_log.py index f492fb18c..e52f93071 100644 --- a/tests/test_gaussian_log.py +++ b/tests/test_gaussian_log.py @@ -74,5 +74,23 @@ def test_forces(self) : def test_virials(self) : self.assertFalse('virials' in self.system.data) + +class TestGaussianLoadPBCLog(unittest.TestCase, TestGaussianLog): + """PBC.""" + def setUp (self) : + self.system = dpdata.LabeledSystem('gaussian/h2pbc.gaussianlog', + fmt = 'gaussian/log') + self.atom_names = ['H'] + self.atom_numbs = [2] + self.nframes = 1 + self.atom_types = [0, 0] + self.cells = (np.eye(3) * 10.0).reshape(1, 3, 3) + + def test_cells(self) : + self.assertTrue(np.allclose(self.system.data['cells'], self.cells)) + + def test_nopbc(self): + self.assertEqual(self.system.nopbc, False) + if __name__ == '__main__': unittest.main()