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