Coverage for /builds/ase/ase/ase/calculators/amber.py: 33.79%
145 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3"""This module defines an ASE interface to Amber16.
5Usage: (Tested only with Amber16, http://ambermd.org/)
7Before usage, input files (infile, topologyfile, incoordfile)
9"""
11import subprocess
13import numpy as np
14from scipy.io import netcdf_file
16import ase.units as units
17from ase.calculators.calculator import Calculator, FileIOCalculator
18from ase.io.amber import read_amber_coordinates, write_amber_coordinates
21class Amber(FileIOCalculator):
22 """Class for doing Amber classical MM calculations.
24 Example:
26 mm.in::
28 Minimization with Cartesian restraints
29 &cntrl
30 imin=1, maxcyc=200, (invoke minimization)
31 ntpr=5, (print frequency)
32 &end
33 """
35 implemented_properties = ['energy', 'forces']
36 discard_results_on_any_change = True
38 def __init__(self, restart=None,
39 ignore_bad_restart_file=FileIOCalculator._deprecated,
40 label='amber', atoms=None, command=None,
41 amber_exe='sander -O ',
42 infile='mm.in', outfile='mm.out',
43 topologyfile='mm.top', incoordfile='mm.crd',
44 outcoordfile='mm_dummy.crd', mdcoordfile=None,
45 **kwargs):
46 """Construct Amber-calculator object.
48 Parameters
49 ==========
50 label: str
51 Name used for all files. May contain a directory.
52 atoms: Atoms object
53 Optional Atoms object to which the calculator will be
54 attached. When restarting, atoms will get its positions and
55 unit-cell updated from file.
56 label: str
57 Prefix to use for filenames (label.in, label.txt, ...).
58 amber_exe: str
59 Name of the amber executable, one can add options like -O
60 and other parameters here
61 infile: str
62 Input filename for amber, contains instuctions about the run
63 outfile: str
64 Logfilename for amber
65 topologyfile: str
66 Name of the amber topology file
67 incoordfile: str
68 Name of the file containing the input coordinates of atoms
69 outcoordfile: str
70 Name of the file containing the output coordinates of atoms
71 this file is not used in case minisation/dynamics is done by ase.
72 It is only relevant
73 if you run MD/optimisation many steps with amber.
75 """
77 self.out = 'mm.log'
79 self.positions = None
80 self.atoms = None
82 self.set(**kwargs)
84 self.amber_exe = amber_exe
85 self.infile = infile
86 self.outfile = outfile
87 self.topologyfile = topologyfile
88 self.incoordfile = incoordfile
89 self.outcoordfile = outcoordfile
90 self.mdcoordfile = mdcoordfile
92 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file,
93 label, atoms, command=command,
94 **kwargs)
96 @property
97 def _legacy_default_command(self):
98 command = (self.amber_exe +
99 ' -i ' + self.infile +
100 ' -o ' + self.outfile +
101 ' -p ' + self.topologyfile +
102 ' -c ' + self.incoordfile +
103 ' -r ' + self.outcoordfile)
104 if self.mdcoordfile is not None:
105 command = command + ' -x ' + self.mdcoordfile
106 return command
108 def write_input(self, atoms=None, properties=None, system_changes=None):
109 """Write updated coordinates to a file."""
111 FileIOCalculator.write_input(self, atoms, properties, system_changes)
112 self.write_coordinates(atoms, self.incoordfile)
114 def read_results(self):
115 """ read energy and forces """
116 self.read_energy()
117 self.read_forces()
119 def write_coordinates(self, atoms, filename):
120 """ write amber coordinates in netCDF format,
121 only rectangular unit cells are allowed"""
122 write_amber_coordinates(atoms, filename)
124 def read_coordinates(self, atoms):
125 """Import AMBER16 netCDF restart files.
127 Reads atom positions and
128 velocities (if available),
129 and unit cell (if available)
131 This may be useful if you have run amber many steps and
132 want to read new positions and velocities
133 """
134 # For historical reasons we edit the input atoms rather than
135 # returning new atoms.
136 _atoms = read_amber_coordinates(self.outcoordfile)
137 atoms.cell[:] = _atoms.cell
138 atoms.pbc[:] = _atoms.pbc
139 atoms.positions[:] = _atoms.positions
140 atoms.set_momenta(_atoms.get_momenta())
142 def read_energy(self, filename='mden'):
143 """ read total energy from amber file """
144 with open(filename, 'r') as fd:
145 lines = fd.readlines()
146 blocks = []
147 while 'L0' in lines[0].split()[0]:
148 blocks.append(lines[0:10])
149 lines = lines[10:]
150 if lines == []:
151 break
152 self.results['energy'] = \
153 float(blocks[-1][6].split()[2]) * units.kcal / units.mol
155 def read_forces(self, filename='mdfrc'):
156 """ read forces from amber file """
157 fd = netcdf_file(filename, 'r', mmap=False)
158 try:
159 forces = fd.variables['forces']
160 self.results['forces'] = forces[-1, :, :] \
161 / units.Ang * units.kcal / units.mol
162 finally:
163 fd.close()
165 def set_charges(self, selection, charges, parmed_filename=None):
166 """ Modify amber topology charges to contain the updated
167 QM charges, needed in QM/MM.
168 Using amber's parmed program to change charges.
169 """
170 qm_list = list(selection)
171 with open(parmed_filename, 'w') as fout:
172 fout.write('# update the following QM charges \n')
173 for i, charge in zip(qm_list, charges):
174 fout.write('change charge @' + str(i + 1) + ' ' +
175 str(charge) + ' \n')
176 fout.write('# Output the topology file \n')
177 fout.write('outparm ' + self.topologyfile + ' \n')
178 parmed_command = ('parmed -O -i ' + parmed_filename +
179 ' -p ' + self.topologyfile +
180 ' > ' + self.topologyfile + '.log 2>&1')
181 subprocess.check_call(parmed_command, shell=True, cwd=self.directory)
183 def get_virtual_charges(self, atoms):
184 with open(self.topologyfile, 'r') as fd:
185 topology = fd.readlines()
186 for n, line in enumerate(topology):
187 if '%FLAG CHARGE' in line:
188 chargestart = n + 2
189 lines1 = topology[chargestart:(chargestart
190 + (len(atoms) - 1) // 5 + 1)]
191 mm_charges = []
192 for line in lines1:
193 for el in line.split():
194 mm_charges.append(float(el) / 18.2223)
195 charges = np.array(mm_charges)
196 return charges
198 def add_virtual_sites(self, positions):
199 return positions # no virtual sites
201 def redistribute_forces(self, forces):
202 return forces
205def map(atoms, top):
206 p = np.zeros((2, len(atoms)), dtype="int")
208 elements = atoms.get_chemical_symbols()
209 unique_elements = np.unique(atoms.get_chemical_symbols())
211 for i in range(len(unique_elements)):
212 idx = 0
213 for j in range(len(atoms)):
214 if elements[j] == unique_elements[i]:
215 idx += 1
216 symbol = unique_elements[i] + np.str(idx)
217 for k in range(len(atoms)):
218 if top.atoms[k].name == symbol:
219 p[0, k] = j
220 p[1, j] = k
221 break
222 return p
225try:
226 import sander
227 have_sander = True
228except ImportError:
229 have_sander = False
232class SANDER(Calculator):
233 """
234 Interface to SANDER using Python interface
236 Requires sander Python bindings from http://ambermd.org/
237 """
238 implemented_properties = ['energy', 'forces']
240 def __init__(self, atoms=None, label=None, top=None, crd=None,
241 mm_options=None, qm_options=None, permutation=None, **kwargs):
242 if not have_sander:
243 raise RuntimeError("sander Python module could not be imported!")
244 Calculator.__init__(self, label, atoms)
245 self.permutation = permutation
246 if qm_options is not None:
247 sander.setup(top, crd.coordinates, crd.box, mm_options, qm_options)
248 else:
249 sander.setup(top, crd.coordinates, crd.box, mm_options)
251 def calculate(self, atoms, properties, system_changes):
252 Calculator.calculate(self, atoms, properties, system_changes)
253 if system_changes:
254 if 'energy' in self.results:
255 del self.results['energy']
256 if 'forces' in self.results:
257 del self.results['forces']
258 if 'energy' not in self.results:
259 if self.permutation is None:
260 crd = np.reshape(atoms.get_positions(), (1, len(atoms), 3))
261 else:
262 crd = np.reshape(atoms.get_positions()
263 [self.permutation[0, :]], (1, len(atoms), 3))
264 sander.set_positions(crd)
265 e, f = sander.energy_forces()
266 self.results['energy'] = e.tot * units.kcal / units.mol
267 if self.permutation is None:
268 self.results['forces'] = (np.reshape(np.array(f),
269 (len(atoms), 3)) *
270 units.kcal / units.mol)
271 else:
272 ff = np.reshape(np.array(f), (len(atoms), 3)) * \
273 units.kcal / units.mol
274 self.results['forces'] = ff[self.permutation[1, :]]
275 if 'forces' not in self.results:
276 if self.permutation is None:
277 crd = np.reshape(atoms.get_positions(), (1, len(atoms), 3))
278 else:
279 crd = np.reshape(atoms.get_positions()[self.permutation[0, :]],
280 (1, len(atoms), 3))
281 sander.set_positions(crd)
282 e, f = sander.energy_forces()
283 self.results['energy'] = e.tot * units.kcal / units.mol
284 if self.permutation is None:
285 self.results['forces'] = (np.reshape(np.array(f),
286 (len(atoms), 3)) *
287 units.kcal / units.mol)
288 else:
289 ff = np.reshape(np.array(f), (len(atoms), 3)) * \
290 units.kcal / units.mol
291 self.results['forces'] = ff[self.permutation[1, :]]