Coverage for /builds/kinetik161/ase/ase/calculators/amber.py: 43.72%

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

14import ase.units as units

15from ase.calculators.calculator import Calculator, FileIOCalculator

18class Amber(FileIOCalculator):

19 """Class for doing Amber classical MM calculations.

21 Example:

23 mm.in::

25 Minimization with Cartesian restraints

26 &cntrl

27 imin=1, maxcyc=200, (invoke minimization)

28 ntpr=5, (print frequency)

29 &end

30 """

32 implemented_properties = ['energy', 'forces']

33 discard_results_on_any_change = True

35 def __init__(self, restart=None,

36 ignore_bad_restart_file=FileIOCalculator._deprecated,

37 label='amber', atoms=None, command=None,

38 amber_exe='sander -O ',

39 infile='mm.in', outfile='mm.out',

40 topologyfile='mm.top', incoordfile='mm.crd',

41 outcoordfile='mm_dummy.crd', mdcoordfile=None,

42 **kwargs):

43 """Construct Amber-calculator object.

45 Parameters

46 ==========

47 label: str

48 Name used for all files. May contain a directory.

49 atoms: Atoms object

50 Optional Atoms object to which the calculator will be

51 attached. When restarting, atoms will get its positions and

52 unit-cell updated from file.

53 label: str

54 Prefix to use for filenames (label.in, label.txt, ...).

55 amber_exe: str

56 Name of the amber executable, one can add options like -O

57 and other parameters here

58 infile: str

59 Input filename for amber, contains instuctions about the run

60 outfile: str

61 Logfilename for amber

62 topologyfile: str

63 Name of the amber topology file

64 incoordfile: str

65 Name of the file containing the input coordinates of atoms

66 outcoordfile: str

67 Name of the file containing the output coordinates of atoms

68 this file is not used in case minisation/dynamics is done by ase.

69 It is only relevant

70 if you run MD/optimisation many steps with amber.

72 """

74 self.out = 'mm.log'

76 self.positions = None

77 self.atoms = None

79 self.set(**kwargs)

81 self.amber_exe = amber_exe

82 self.infile = infile

83 self.outfile = outfile

84 self.topologyfile = topologyfile

85 self.incoordfile = incoordfile

86 self.outcoordfile = outcoordfile

87 self.mdcoordfile = mdcoordfile

89 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file,

90 label, atoms, command=command,

91 **kwargs)

93 @property

94 def _legacy_default_command(self):

95 command = (self.amber_exe +

96 ' -i ' + self.infile +

97 ' -o ' + self.outfile +

98 ' -p ' + self.topologyfile +

99 ' -c ' + self.incoordfile +

100 ' -r ' + self.outcoordfile)

101 if self.mdcoordfile is not None:

102 command = command + ' -x ' + self.mdcoordfile

103 return command

104

105 def write_input(self, atoms=None, properties=None, system_changes=None):

106 """Write updated coordinates to a file."""

107

108 FileIOCalculator.write_input(self, atoms, properties, system_changes)

109 self.write_coordinates(atoms)

110

111 def read_results(self):

112 """ read energy and forces """

113 self.read_energy()

114 self.read_forces()

115

116 def write_coordinates(self, atoms, filename=''):

117 """ write amber coordinates in netCDF format,

118 only rectangular unit cells are allowed"""

119 if filename == '':

120 filename = self.incoordfile

121 from scipy.io import netcdf_file

122 fout = netcdf_file(filename, 'w')

123 # dimension

124 fout.Conventions = 'AMBERRESTART'

125 fout.ConventionVersion = "1.0"

126 fout.title = 'Ase-generated-amber-restart-file'

127 fout.application = "AMBER"

128 fout.program = "ASE"

129 fout.programVersion = "1.0"

130 fout.createDimension('cell_spatial', 3)

131 fout.createDimension('label', 5)

132 fout.createDimension('cell_angular', 3)

133 fout.createDimension('time', 1)

134 time = fout.createVariable('time', 'd', ('time',))

135 time.units = 'picosecond'

136 time[0] = 0

137 fout.createDimension('spatial', 3)

138 spatial = fout.createVariable('spatial', 'c', ('spatial',))

139 spatial[:] = np.asarray(list('xyz'))

140 # spatial = 'xyz'

141

142 natom = len(atoms)

143 fout.createDimension('atom', natom)

144 coordinates = fout.createVariable('coordinates', 'd',

145 ('atom', 'spatial'))

146 coordinates.units = 'angstrom'

147 coordinates[:] = atoms.get_positions()[:]

148

149 if atoms.get_velocities() is not None:

150 velocities = fout.createVariable('velocities', 'd',

151 ('atom', 'spatial'))

152 velocities.units = 'angstrom/picosecond'

153 velocities[:] = atoms.get_velocities()[:]

154

155 # title

156 cell_angular = fout.createVariable('cell_angular', 'c',

157 ('cell_angular', 'label'))

158 cell_angular[0] = np.asarray(list('alpha'))

159 cell_angular[1] = np.asarray(list('beta '))

160 cell_angular[2] = np.asarray(list('gamma'))

161

162 # title

163 cell_spatial = fout.createVariable('cell_spatial', 'c',

164 ('cell_spatial',))

165 cell_spatial[0], cell_spatial[1], cell_spatial[2] = 'a', 'b', 'c'

166

167 # data

168 cell_lengths = fout.createVariable('cell_lengths', 'd',

169 ('cell_spatial',))

170 cell_lengths.units = 'angstrom'

171 cell_lengths[0] = atoms.get_cell()[0, 0]

172 cell_lengths[1] = atoms.get_cell()[1, 1]

173 cell_lengths[2] = atoms.get_cell()[2, 2]

174

175 cell_angles = fout.createVariable('cell_angles', 'd',

176 ('cell_angular',))

177 box_alpha, box_beta, box_gamma = 90.0, 90.0, 90.0

178 cell_angles[0] = box_alpha

179 cell_angles[1] = box_beta

180 cell_angles[2] = box_gamma

181

182 cell_angles.units = 'degree'

183 fout.close()

184

185 def read_coordinates(self, atoms, filename=''):

186 """Import AMBER16 netCDF restart files.

187

188 Reads atom positions and

189 velocities (if available),

190 and unit cell (if available)

191

192 This may be useful if you have run amber many steps and

193 want to read new positions and velocities

194 """

195

196 if filename == '':

197 filename = self.outcoordfile

198

199 import numpy as np

200 from scipy.io import netcdf

201

202 import ase.units as units

203

204 fin = netcdf.netcdf_file(filename, 'r')

205 all_coordinates = fin.variables['coordinates'][:]

206 get_last_frame = False

207 if hasattr(all_coordinates, 'ndim'):

208 if all_coordinates.ndim == 3:

209 get_last_frame = True

210 elif hasattr(all_coordinates, 'shape'):

211 if len(all_coordinates.shape) == 3:

212 get_last_frame = True

213 if get_last_frame:

214 all_coordinates = all_coordinates[-1]

215 atoms.set_positions(all_coordinates)

216 if 'velocities' in fin.variables:

217 all_velocities = fin.variables['velocities'][:] / (1000 * units.fs)

218 if get_last_frame:

219 all_velocities = all_velocities[-1]

220 atoms.set_velocities(all_velocities)

221 if 'cell_lengths' in fin.variables:

222 all_abc = fin.variables['cell_lengths']

223 if get_last_frame:

224 all_abc = all_abc[-1]

225 a, b, c = all_abc

226 all_angles = fin.variables['cell_angles']

227 if get_last_frame:

228 all_angles = all_angles[-1]

229 alpha, beta, gamma = all_angles

230

231 if (all(angle > 89.99 for angle in [alpha, beta, gamma]) and

232 all(angle < 90.01 for angle in [alpha, beta, gamma])):

233 atoms.set_cell(

234 np.array([[a, 0, 0],

235 [0, b, 0],

236 [0, 0, c]]))

237 atoms.set_pbc(True)

238 else:

239 raise NotImplementedError('only rectangular cells are'

240 ' implemented in ASE-AMBER')

241

242 else:

243 atoms.set_pbc(False)

244

245 def read_energy(self, filename='mden'):

246 """ read total energy from amber file """

247 with open(filename) as fd:

248 lines = fd.readlines()

249 self.results['energy'] = \

250 float(lines[16].split()[2]) * units.kcal / units.mol

251

252 def read_forces(self, filename='mdfrc'):

253 """ read forces from amber file """

254 fd = netcdf.netcdf_file(filename, 'r')

255 try:

256 forces = fd.variables['forces']

257 self.results['forces'] = forces[-1, :, :] \

258 / units.Ang * units.kcal / units.mol

259 finally:

260 fd.close()

261

262 def set_charges(self, selection, charges, parmed_filename=None):

263 """ Modify amber topology charges to contain the updated

264 QM charges, needed in QM/MM.

265 Using amber's parmed program to change charges.

266 """

267 qm_list = list(selection)

268 with open(parmed_filename, 'w') as fout:

269 fout.write('# update the following QM charges \n')

270 for i, charge in zip(qm_list, charges):

271 fout.write('change charge @' + str(i + 1) + ' ' +

272 str(charge) + ' \n')

273 fout.write('# Output the topology file \n')

274 fout.write('outparm ' + self.topologyfile + ' \n')

275 parmed_command = ('parmed -O -i ' + parmed_filename +

276 ' -p ' + self.topologyfile +

277 ' > ' + self.topologyfile + '.log 2>&1')

278 subprocess.check_call(parmed_command, shell=True, cwd=self.directory)

279

280 def get_virtual_charges(self, atoms):

281 with open(self.topologyfile) as fd:

282 topology = fd.readlines()

283 for n, line in enumerate(topology):

284 if '%FLAG CHARGE' in line:

285 chargestart = n + 2

286 lines1 = topology[chargestart:(chargestart

287 + (len(atoms) - 1) // 5 + 1)]

288 mm_charges = []

289 for line in lines1:

290 for el in line.split():

291 mm_charges.append(float(el) / 18.2223)

292 charges = np.array(mm_charges)

293 return charges

294

295 def add_virtual_sites(self, positions):

296 return positions # no virtual sites

297

298 def redistribute_forces(self, forces):

299 return forces

300

301

302def map(atoms, top):

303 p = np.zeros((2, len(atoms)), dtype="int")

304

305 elements = atoms.get_chemical_symbols()

306 unique_elements = np.unique(atoms.get_chemical_symbols())

307

308 for i in range(len(unique_elements)):

309 idx = 0

310 for j in range(len(atoms)):

311 if elements[j] == unique_elements[i]:

312 idx += 1

313 symbol = unique_elements[i] + np.str(idx)

314 for k in range(len(atoms)):

315 if top.atoms[k].name == symbol:

316 p[0, k] = j

317 p[1, j] = k

318 break

319 return p

320

321

322try:

323 import sander

324 have_sander = True

325except ImportError:

326 have_sander = False

327

328

329class SANDER(Calculator):

330 """

331 Interface to SANDER using Python interface

332

333 Requires sander Python bindings from http://ambermd.org/

334 """

335 implemented_properties = ['energy', 'forces']

336

337 def __init__(self, atoms=None, label=None, top=None, crd=None,

338 mm_options=None, qm_options=None, permutation=None, **kwargs):

339 if not have_sander:

340 raise RuntimeError("sander Python module could not be imported!")

341 Calculator.__init__(self, label, atoms)

342 self.permutation = permutation

343 if qm_options is not None:

344 sander.setup(top, crd.coordinates, crd.box, mm_options, qm_options)

345 else:

346 sander.setup(top, crd.coordinates, crd.box, mm_options)

347

348 def calculate(self, atoms, properties, system_changes):

349 Calculator.calculate(self, atoms, properties, system_changes)

350 if system_changes:

351 if 'energy' in self.results:

352 del self.results['energy']

353 if 'forces' in self.results:

354 del self.results['forces']

355 if 'energy' not in self.results:

356 if self.permutation is None:

357 crd = np.reshape(atoms.get_positions(), (1, len(atoms), 3))

358 else:

359 crd = np.reshape(atoms.get_positions()

360 [self.permutation[0, :]], (1, len(atoms), 3))

361 sander.set_positions(crd)

362 e, f = sander.energy_forces()

363 self.results['energy'] = e.tot * units.kcal / units.mol

364 if self.permutation is None:

365 self.results['forces'] = (np.reshape(np.array(f),

366 (len(atoms), 3)) *

367 units.kcal / units.mol)

368 else:

369 ff = np.reshape(np.array(f), (len(atoms), 3)) * \

370 units.kcal / units.mol

371 self.results['forces'] = ff[self.permutation[1, :]]

372 if 'forces' not in self.results:

373 if self.permutation is None:

374 crd = np.reshape(atoms.get_positions(), (1, len(atoms), 3))

375 else:

376 crd = np.reshape(atoms.get_positions()[self.permutation[0, :]],

377 (1, len(atoms), 3))

378 sander.set_positions(crd)

379 e, f = sander.energy_forces()

380 self.results['energy'] = e.tot * units.kcal / units.mol

381 if self.permutation is None:

382 self.results['forces'] = (np.reshape(np.array(f),

383 (len(atoms), 3)) *

384 units.kcal / units.mol)

385 else:

386 ff = np.reshape(np.array(f), (len(atoms), 3)) * \

387 units.kcal / units.mol

388 self.results['forces'] = ff[self.permutation[1, :]]