Coverage for cclib/bridge/cclib2ase.py : 95%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# -*- coding: utf-8 -*-
2#
3# Copyright (c) 2020, the cclib development team
4#
5# This file is part of cclib (http://cclib.github.io) and is distributed under
6# the terms of the BSD 3-Clause License.
8"""Bridge for using cclib data in ASE (https://wiki.fysik.dtu.dk/ase/)."""
10import numpy as np
12from cclib.parser.data import ccData
13from cclib.parser.utils import find_package
15_found_ase = find_package("ase")
16if _found_ase:
17 from ase import Atoms, units
18 from ase.io.trajectory import Trajectory
19 from ase.calculators.calculator import PropertyNotImplementedError
22def _check_ase(found_ase):
23 if not found_ase:
24 raise ImportError("You must install `ase` to use this function")
27def makease(
28 atomcoords, atomnos, atomcharges=None, atomspins=None, atommasses=None
29):
30 """Create an ASE Atoms object from cclib attributes.
32 ASE requires atomic partial charges and atomic spin densities rather than
33 molecular charge and multiplicity, so we follow how other interfaces have
34 done (e.g., MOPAC, Gaussian and XTB) and require atomcharges and atomspins,
35 or leave undefined.
37 Inputs:
38 atomcoords - two dimensional array-like with atomic coordinates.
39 atomnos - one dimensional array-like with atomic numbers.
40 atomcharges - one dimensional array-like with atomic charges.
41 atomspins - one dimensional array-like with atomic spin densities.
42 atommasses - one dimensional array-like with atomic masses.
43 """
44 _check_ase(_found_ase)
45 return Atoms(
46 positions=atomcoords,
47 numbers=atomnos,
48 masses=atommasses,
49 charges=atomcharges,
50 magmoms=atomspins,
51 )
54def write_trajectory(filename, ccdata, popname="mulliken", index=None):
55 """Write an ASE Trajectory object from a ccData object.
57 We try to write the following properties: atomcoords, atomnos, atomcharges,
58 atomspins, atommasses, scfenergies, grads, moments (dipole only) and
59 freeenergy. No charge or mult is written since ASE calculates them from
60 atomcharges and atomspins.
62 When a program do a final single point calculation at the end of an
63 optimization (e.g., ORCA), we'll have unequal number of grads and
64 atomcoords. This means the last geometry in the written .traj will lack
65 forces.
67 Some unit conversions are done here, but things are converted back in
68 read_trajectory/makecclib.
70 Inputs:
71 filename - path to traj file to be written.
72 ccdata - an instance of ccData.
73 popname - population analysis to use for atomic partial charges and
74 atomic spin densities. Molecular charge and multiplicity are
75 evaluated from them.
76 index - sequence of integers indicating which atomcoords indices should
77 be exported. By default, all are exported.
78 """
79 _check_ase(_found_ase)
80 traj = Trajectory(filename, "w")
81 for i, atomcoords in enumerate(ccdata.atomcoords):
82 if index is not None and i not in index:
83 continue
85 atomspins = None
86 if hasattr(ccdata, "atomspins"):
87 atomspins = ccdata.atomspins[popname]
88 atoms = makease(
89 atomcoords,
90 ccdata.atomnos,
91 ccdata.atomcharges[popname],
92 atomspins,
93 ccdata.atommasses,
94 )
96 properties = {}
97 if hasattr(ccdata, "scfenergies"):
98 properties.update({"energy": ccdata.scfenergies[i]})
99 if hasattr(ccdata, "grads"):
100 try:
101 properties.update(
102 {"forces": -ccdata.grads[i] * units.Hartree / units.Bohr}
103 )
104 except IndexError:
105 pass
107 if i == len(ccdata.atomcoords) - 1: # last geometry
108 if hasattr(ccdata, "moments"):
109 properties.update({"dipole": ccdata.moments[1] * units.Bohr})
110 if hasattr(ccdata, "free_energy"):
111 properties.update(
112 {"free_energy": ccdata.freeenergy * units.Hartree}
113 )
115 traj.write(atoms, **properties)
118def read_trajectory(filename):
119 """Read an ASE Trajectory object and return a ccData object.
121 The returned object has everything write_trajectory writes, plus natom,
122 charge, mult and temperature.
124 The following properties are taken from the last frame: atomnos,
125 atomcharges, atomspins, atommasses, moments, freeenergy and temperature.
126 charge, mult and natom also represent the last frame, since they depend on
127 other propertes read from the last frame.
129 Bear in mind that ASE calculates temperature from the kinetic energy, so
130 anything "static" (which includes anything cclib parses) will have zero
131 temperature.
133 Inputs:
134 filename - path to traj file to be read.
135 """
136 _check_ase(_found_ase)
137 attributes = {"atomcoords": [], "scfenergies": [], "grads": []}
139 for atoms in Trajectory(filename, "r"):
140 ccdata = makecclib(atoms)
141 attributes["atomcoords"].append(ccdata.atomcoords[-1])
142 if hasattr(ccdata, "scfenergies"):
143 attributes["scfenergies"].append(ccdata.scfenergies[-1])
144 if hasattr(ccdata, "grads"):
145 attributes["grads"].append(ccdata.grads[-1])
147 # ccdata is now last frame
148 attributes["atomnos"] = ccdata.atomnos
149 attributes["atomcharges"] = ccdata.atomcharges
150 attributes["atomspins"] = ccdata.atomspins
151 attributes["atommasses"] = ccdata.atommasses
153 if hasattr(ccdata, "moments"):
154 attributes["moments"] = ccdata.moments
155 if hasattr(ccdata, "freeenergy"):
156 attributes["freeenergy"] = ccdata.freeenergy
158 # remove if empty
159 if not attributes["scfenergies"]:
160 del attributes["scfenergies"]
161 if not attributes["grads"]:
162 del attributes["grads"]
164 # extra stuff we can't write in write_trajectory
165 attributes["temperature"] = ccdata.temperature
166 attributes["charge"] = ccdata.charge
167 attributes["mult"] = ccdata.mult
168 attributes["natom"] = ccdata.natom
170 return ccData(attributes)
173def makecclib(atoms, popname="mulliken"):
174 """Create cclib attributes and return a ccData from an ASE Atoms object.
176 Available data (such as forces/gradients and potential energy/free
177 energy) is assumed to be from SCF (see
178 https://wiki.fysik.dtu.dk/ase/ase/atoms.html#adding-a-calculator).
180 Bear in mind that ASE calculates temperature from the kinetic energy, so
181 anything "static" (which includes anything cclib parses) will have zero
182 temperature.
184 Inputs:
185 atoms - an instance of ASE `Atoms`
186 popname - population analysis to use for atomic partial charges and
187 atomic spin densities. Molecular charge and multiplicity are
188 evaluated from them.
189 """
190 _check_ase(_found_ase)
191 attributes = {
192 "atomcoords": np.array([atoms.get_positions()]),
193 "atomnos": atoms.get_atomic_numbers(),
194 "atommasses": atoms.get_masses(),
195 "natom": atoms.get_global_number_of_atoms(),
196 }
198 try:
199 attributes["atomcharges"] = {popname: atoms.get_charges()}
200 except (PropertyNotImplementedError, RuntimeError):
201 attributes["atomcharges"] = {popname: atoms.get_initial_charges()}
202 try:
203 attributes["atomspins"] = {popname: atoms.get_magnetic_moments()}
204 except (PropertyNotImplementedError, RuntimeError):
205 attributes["atomspins"] = {
206 popname: atoms.get_initial_magnetic_moments()
207 }
209 # the following is how ASE determines charge and multiplicity from initial
210 # charges and initial magnetic moments in its Gaussian interface
211 # (https://gitlab.com/ase/ase/-/blob/a26bda2160527ca7afc0135c69e4367a5bc5a264/ase/io/gaussian.py#L105)
212 attributes["charge"] = attributes["atomcharges"][popname].sum()
213 attributes["mult"] = attributes["atomspins"][popname].sum() + 1
215 try:
216 attributes["scfenergies"] = np.array([atoms.get_potential_energy()])
217 except RuntimeError:
218 pass
219 try:
220 attributes["grads"] = (
221 -np.array([atoms.get_forces()]) * units.Bohr / units.Hartree
222 )
223 except RuntimeError:
224 pass
225 try:
226 attributes["moments"] = [
227 atoms.get_center_of_mass(),
228 atoms.get_dipole_moment() / units.Bohr,
229 ]
230 except RuntimeError:
231 pass
232 try:
233 attributes["freeenergy"] = (
234 atoms.get_potential_energy(force_consistent=True) / units.Hartree
235 )
236 except RuntimeError:
237 pass
239 attributes["temperature"] = atoms.get_temperature()
241 return ccData(attributes)
244del find_package