forked from lijiext/lammps
727 lines
27 KiB
Python
727 lines
27 KiB
Python
|
#!/usr/bin/env python
|
||
|
|
# polarizer.py - add Drude oscillators to LAMMPS data file.
|
||
|
|
# Agilio Padua <agilio.padua@univ-bpclermont.fr>
|
||
|
|
# Alain Dequidt <alain.dequidt@univ-bpclermont.fr>
|
||
|
|
# version 2015/07/17
|
||
|
|
|
||
|
|
import sys
|
||
|
|
import argparse
|
||
|
|
import random
|
||
|
|
from copy import deepcopy
|
||
|
|
|
||
|
|
usage = """Add Drude oscillators to LAMMPS data file.
|
||
|
|
|
||
|
|
Format of file containing specification of Drude oscillators:
|
||
|
|
|
||
|
|
# type dm/u dq/e k/(kJ/molA2) alpha/A3 thole
|
||
|
|
C3H 1.0 0.0 4184.0 2.051 2.6
|
||
|
|
...
|
||
|
|
|
||
|
|
* dm is the mass to place on the Drude particle (taken from its core),
|
||
|
|
* dq is the charge to place on the Drude particle (taken from its core),
|
||
|
|
* k is the harmonic force constant of the bond between core and Drude,
|
||
|
|
* alpha is the polarizability,
|
||
|
|
* thole is a parameter of the Thole damping function.
|
||
|
|
|
||
|
|
A Drude particle is created for each atom in the LAMMPS data file
|
||
|
|
that corresponds to an atom type given in the Drude file.
|
||
|
|
Since LAMMPS uses numbers for atom types in the data file, a comment
|
||
|
|
after each line in the Masses section has to be introduced to allow
|
||
|
|
identification of the atom types within the force field database:
|
||
|
|
|
||
|
|
Masses
|
||
|
|
|
||
|
|
1 12.011 # C3H
|
||
|
|
2 12.011 # CTO
|
||
|
|
...
|
||
|
|
|
||
|
|
This script will add new atom types, new bond types, new atoms and
|
||
|
|
new bonds to the data file.
|
||
|
|
|
||
|
|
It will also print some commands to be included in the LAMMPS input script,
|
||
|
|
which are related to the topology and force field, namely fix drude,
|
||
|
|
pair_style and pair coeff_commands. For information on thermostating please
|
||
|
|
read the documentation of the USER-DRUDE package.
|
||
|
|
|
||
|
|
This tool can also be used to revert a Drude-polarized data file to a
|
||
|
|
non-polarizable one.
|
||
|
|
"""
|
||
|
|
|
||
|
|
# keywords of header and main sections (from data.py in Pizza.py)
|
||
|
|
|
||
|
|
hkeywords = ["atoms", "ellipsoids", "lines", "triangles", "bodies",
|
||
|
|
"bonds", "angles", "dihedrals", "impropers",
|
||
|
|
"atom types", "bond types", "angle types", "dihedral types",
|
||
|
|
"improper types", "xlo xhi", "ylo yhi", "zlo zhi", "xy xz yz"]
|
||
|
|
|
||
|
|
skeywords = [["Masses", "atom types"],
|
||
|
|
["Pair Coeffs", "atom types"],
|
||
|
|
["Bond Coeffs", "bond types"], ["Angle Coeffs", "angle types"],
|
||
|
|
["Dihedral Coeffs", "dihedral types"],
|
||
|
|
["Improper Coeffs", "improper types"],
|
||
|
|
["BondBond Coeffs", "angle types"],
|
||
|
|
["BondAngle Coeffs", "angle types"],
|
||
|
|
["MiddleBondTorsion Coeffs", "dihedral types"],
|
||
|
|
["EndBondTorsion Coeffs", "dihedral types"],
|
||
|
|
["AngleTorsion Coeffs", "dihedral types"],
|
||
|
|
["AngleAngleTorsion Coeffs", "dihedral types"],
|
||
|
|
["BondBond13 Coeffs", "dihedral types"],
|
||
|
|
["AngleAngle Coeffs", "improper types"],
|
||
|
|
["Atoms", "atoms"], ["Velocities", "atoms"],
|
||
|
|
["Ellipsoids", "ellipsoids"],
|
||
|
|
["Lines", "lines"], ["Triangles", "triangles"],
|
||
|
|
["Bodies", "bodies"],
|
||
|
|
["Bonds", "bonds"],
|
||
|
|
["Angles", "angles"], ["Dihedrals", "dihedrals"],
|
||
|
|
["Impropers", "impropers"], ["Molecules", "atoms"]]
|
||
|
|
|
||
|
|
|
||
|
|
def massline(att):
|
||
|
|
return "{0:4d} {1:8.3f} # {2}\n".format(att['id'], att['m'], att['type'])
|
||
|
|
|
||
|
|
def bdtline(bdt):
|
||
|
|
return "{0:4d} {1:12.6f} {2:12.6f} {3}\n".format(bdt['id'], bdt['k'],
|
||
|
|
bdt['r0'], bdt['note'])
|
||
|
|
|
||
|
|
def atomline(at):
|
||
|
|
return "{0:7d} {1:7d} {2:4d} {3:8.4f} {4:13.6e} {5:13.6e} {6:13.6e} "\
|
||
|
|
" {7}\n".format(at['n'], at['mol'], at['id'], at['q'],
|
||
|
|
at['x'], at['y'], at['z'], at['note'])
|
||
|
|
|
||
|
|
def bondline(bd):
|
||
|
|
return "{0:7d} {1:4d} {2:7d} {3:7d} {4}\n".format(bd['n'], bd['id'],
|
||
|
|
bd['i'], bd['j'], bd['note'])
|
||
|
|
|
||
|
|
def velline(at):
|
||
|
|
return "{0:7d} {1:13.6e} {2:13.6e} {3:13.6e} \n".format(at['n'],
|
||
|
|
at['vx'], at['vy'], at['vz'])
|
||
|
|
|
||
|
|
# --------------------------------------
|
||
|
|
|
||
|
|
|
||
|
|
class Data(object):
|
||
|
|
|
||
|
|
def __init__(self, datafile):
|
||
|
|
'''read LAMMPS data file (from data.py in Pizza.py)'''
|
||
|
|
|
||
|
|
# for extract method
|
||
|
|
self.atomtypes = []
|
||
|
|
self.bondtypes = []
|
||
|
|
self.atoms = []
|
||
|
|
self.bonds = []
|
||
|
|
self.idmap = {}
|
||
|
|
|
||
|
|
self.nselect = 1
|
||
|
|
|
||
|
|
f = open(datafile, "r")
|
||
|
|
|
||
|
|
self.title = f.readline()
|
||
|
|
self.names = {}
|
||
|
|
|
||
|
|
headers = {}
|
||
|
|
while 1:
|
||
|
|
line = f.readline().strip()
|
||
|
|
if '#' in line:
|
||
|
|
line = line[:line.index('#')].strip()
|
||
|
|
if len(line) == 0:
|
||
|
|
continue
|
||
|
|
found = 0
|
||
|
|
for keyword in hkeywords:
|
||
|
|
if keyword in line:
|
||
|
|
found = 1
|
||
|
|
words = line.split()
|
||
|
|
if keyword == "xlo xhi" or keyword == "ylo yhi" or \
|
||
|
|
keyword == "zlo zhi":
|
||
|
|
headers[keyword] = (float(words[0]), float(words[1]))
|
||
|
|
elif keyword == "xy xz yz":
|
||
|
|
headers[keyword] = \
|
||
|
|
(float(words[0]), float(words[1]), float(words[2]))
|
||
|
|
else:
|
||
|
|
headers[keyword] = int(words[0])
|
||
|
|
if not found:
|
||
|
|
break
|
||
|
|
|
||
|
|
sections = {}
|
||
|
|
while 1:
|
||
|
|
if len(line) > 0:
|
||
|
|
found = 0
|
||
|
|
for pair in skeywords:
|
||
|
|
keyword, length = pair[0], pair[1]
|
||
|
|
if keyword == line:
|
||
|
|
found = 1
|
||
|
|
if length not in headers:
|
||
|
|
raise RuntimeError("data section {} "\
|
||
|
|
"has no matching header value".format(line))
|
||
|
|
f.readline()
|
||
|
|
list_ = []
|
||
|
|
for _ in range(headers[length]):
|
||
|
|
list_.append(f.readline())
|
||
|
|
sections[keyword] = list_
|
||
|
|
if not found:
|
||
|
|
raise RuntimeError("invalid section {} in data"\
|
||
|
|
" file".format(line))
|
||
|
|
#f.readline()
|
||
|
|
line = f.readline()
|
||
|
|
if not line:
|
||
|
|
break
|
||
|
|
if '#' in line:
|
||
|
|
line = line[:line.index('#')]
|
||
|
|
line = line.strip()
|
||
|
|
|
||
|
|
f.close()
|
||
|
|
self.headers = headers
|
||
|
|
self.sections = sections
|
||
|
|
|
||
|
|
def write(self, filename):
|
||
|
|
'''write out a LAMMPS data file (from data.py in Pizza.py)'''
|
||
|
|
|
||
|
|
with open(filename, "w") as f:
|
||
|
|
f.write(self.title + '\n')
|
||
|
|
for keyword in hkeywords:
|
||
|
|
if keyword in self.headers:
|
||
|
|
if keyword == "xlo xhi" or keyword == "ylo yhi" or \
|
||
|
|
keyword == "zlo zhi":
|
||
|
|
f.write("{0:f} {1:f} {2}\n".format(
|
||
|
|
self.headers[keyword][0],
|
||
|
|
self.headers[keyword][1], keyword))
|
||
|
|
elif keyword == "xy xz yz":
|
||
|
|
f.write("{0:f} {1:f} {2:f} {3}\n".format(
|
||
|
|
self.headers[keyword][0],
|
||
|
|
self.headers[keyword][1],
|
||
|
|
self.headers[keyword][2], keyword))
|
||
|
|
else:
|
||
|
|
f.write("{0:d} {1}\n".format(self.headers[keyword],
|
||
|
|
keyword))
|
||
|
|
for pair in skeywords:
|
||
|
|
keyword = pair[0]
|
||
|
|
if keyword in self.sections:
|
||
|
|
f.write("\n{}\n\n".format(keyword))
|
||
|
|
for line in self.sections[keyword]:
|
||
|
|
f.write(line)
|
||
|
|
|
||
|
|
def extract_nonpol(self):
|
||
|
|
"""extract atom and bond info from nonpolarizable data"""
|
||
|
|
|
||
|
|
# extract atom IDs
|
||
|
|
missinglabels = False
|
||
|
|
for line in self.sections['Masses']:
|
||
|
|
tok = line.split()
|
||
|
|
if len(tok) < 4:
|
||
|
|
print("error: missing type for atom ID " + tok[0] +
|
||
|
|
" in Masses section")
|
||
|
|
missinglabels = True
|
||
|
|
continue
|
||
|
|
atomtype = {}
|
||
|
|
atomtype['id'] = int(tok[0])
|
||
|
|
atomtype['m'] = float(tok[1])
|
||
|
|
atomtype['type'] = tok[3]
|
||
|
|
self.atomtypes.append(atomtype)
|
||
|
|
|
||
|
|
if missinglabels:
|
||
|
|
sys.exit(0)
|
||
|
|
|
||
|
|
# extract atom registers
|
||
|
|
for line in self.sections['Atoms']:
|
||
|
|
tok = line.split()
|
||
|
|
atom = {}
|
||
|
|
atom['n'] = int(tok[0])
|
||
|
|
atom['mol'] = int(tok[1])
|
||
|
|
atom['id'] = int(tok[2])
|
||
|
|
atom['q'] = float(tok[3])
|
||
|
|
atom['x'] = float(tok[4])
|
||
|
|
atom['y'] = float(tok[5])
|
||
|
|
atom['z'] = float(tok[6])
|
||
|
|
#atom['note'] = ''.join([s + ' ' for s in tok[7:]]).strip()
|
||
|
|
atom['note'] = ' '.join(tok[7:])
|
||
|
|
self.atoms.append(atom)
|
||
|
|
self.idmap[atom['n']] = atom
|
||
|
|
|
||
|
|
if 'Velocities' in self.sections:
|
||
|
|
for line in self.sections['Velocities']:
|
||
|
|
tok = line.split()
|
||
|
|
atom = self.idmap[int(tok[0])]
|
||
|
|
atom['vx'] = float(tok[1])
|
||
|
|
atom['vy'] = float(tok[2])
|
||
|
|
atom['vz'] = float(tok[3])
|
||
|
|
|
||
|
|
def polarize(self, drude):
|
||
|
|
"""add Drude particles"""
|
||
|
|
|
||
|
|
if 'Pair Coeffs' in self.sections:
|
||
|
|
raise RuntimeError("cannot polarize a data with Pair Coeffs")
|
||
|
|
|
||
|
|
self.extract_nonpol()
|
||
|
|
|
||
|
|
natom = self.headers['atoms']
|
||
|
|
nbond = self.headers['bonds']
|
||
|
|
nattype = self.headers['atom types']
|
||
|
|
nbdtype = self.headers['bond types']
|
||
|
|
|
||
|
|
# create new atom types (IDs) for Drude particles and modify cores
|
||
|
|
newattypes = []
|
||
|
|
for att in self.atomtypes:
|
||
|
|
att['dflag'] = 'n'
|
||
|
|
for ddt in drude.types:
|
||
|
|
if ddt['type'] == att['type']:
|
||
|
|
nattype += 1
|
||
|
|
newid = {}
|
||
|
|
newid['id'] = ddt['id'] = nattype
|
||
|
|
newid['m'] = ddt['dm']
|
||
|
|
att['m'] -= ddt['dm']
|
||
|
|
# label drude particles and cores
|
||
|
|
att['dflag'] = 'c'
|
||
|
|
newid['dflag'] = 'd'
|
||
|
|
newid['type'] = att['type'] + ' DP'
|
||
|
|
att['type'] += ' DC'
|
||
|
|
ddt['type'] += ' DC'
|
||
|
|
newattypes.append(newid)
|
||
|
|
break
|
||
|
|
|
||
|
|
self.headers['atom types'] += len(newattypes)
|
||
|
|
self.sections['Masses'] = []
|
||
|
|
for att in self.atomtypes + newattypes:
|
||
|
|
self.sections['Masses'].append(massline(att))
|
||
|
|
|
||
|
|
# create new bond types for core-drude bonds
|
||
|
|
newbdtypes = []
|
||
|
|
for att in self.atomtypes:
|
||
|
|
for ddt in drude.types:
|
||
|
|
if ddt['type'] == att['type']:
|
||
|
|
nbdtype += 1
|
||
|
|
newbdtype = {}
|
||
|
|
newbdtype['id'] = ddt['bdid'] = nbdtype
|
||
|
|
newbdtype['k'] = ddt['k']
|
||
|
|
newbdtype['r0'] = 0.0
|
||
|
|
newbdtype['note'] = '# ' + ddt['type'] + '-DP'
|
||
|
|
newbdtypes.append(newbdtype)
|
||
|
|
break
|
||
|
|
|
||
|
|
self.headers['bond types'] += len(newbdtypes)
|
||
|
|
for bdt in newbdtypes:
|
||
|
|
self.sections['Bond Coeffs'].append(bdtline(bdt))
|
||
|
|
|
||
|
|
# create new atoms for Drude particles and new bonds with their cores
|
||
|
|
random.seed(123)
|
||
|
|
newatoms = []
|
||
|
|
newbonds = []
|
||
|
|
for atom in self.atoms:
|
||
|
|
atom['dflag'] = '' # [c]ore, [d]rude, [n]on-polarizable
|
||
|
|
atom['dd'] = 0 # partner drude or core
|
||
|
|
for att in self.atomtypes:
|
||
|
|
if att['id'] == atom['id']:
|
||
|
|
break
|
||
|
|
for ddt in drude.types:
|
||
|
|
if ddt['type'] == att['type']:
|
||
|
|
natom += 1
|
||
|
|
newatom = deepcopy(atom)
|
||
|
|
newatom['n'] = natom
|
||
|
|
self.idmap[natom] = newatom
|
||
|
|
newatom['id'] = ddt['id']
|
||
|
|
newatom['q'] = ddt['dq']
|
||
|
|
newatom['note'] = atom['note']
|
||
|
|
if '#' not in newatom['note']:
|
||
|
|
newatom['note'] += ' #'
|
||
|
|
newatom['note'] += ' DP'
|
||
|
|
newatom['dflag'] = 'd'
|
||
|
|
newatom['dd'] = atom['n']
|
||
|
|
|
||
|
|
# avoid superposition of cores and Drudes
|
||
|
|
newatom['x'] += 0.1 * (random.random() - 0.5)
|
||
|
|
newatom['y'] += 0.1 * (random.random() - 0.5)
|
||
|
|
newatom['z'] += 0.1 * (random.random() - 0.5)
|
||
|
|
if 'Velocities' in self.sections:
|
||
|
|
newatom['vx'] = atom['vx']
|
||
|
|
newatom['vy'] = atom['vy']
|
||
|
|
newatom['vz'] = atom['vz']
|
||
|
|
|
||
|
|
newatoms.append(newatom)
|
||
|
|
atom['q'] -= ddt['dq']
|
||
|
|
atom['dflag'] = 'c'
|
||
|
|
atom['dd'] = natom
|
||
|
|
if '#' not in atom['note']:
|
||
|
|
atom['note'] += ' #'
|
||
|
|
atom['note'] += ' DC'
|
||
|
|
|
||
|
|
nbond += 1
|
||
|
|
newbond = {}
|
||
|
|
newbond['n'] = nbond
|
||
|
|
newbond['id'] = ddt['bdid']
|
||
|
|
newbond['i'] = atom['n']
|
||
|
|
newbond['j'] = newatom['n']
|
||
|
|
newbond['note'] = '# ' + ddt['type'] + '-DP'
|
||
|
|
newbonds.append(newbond)
|
||
|
|
break
|
||
|
|
|
||
|
|
self.headers['atoms'] += len(newatoms)
|
||
|
|
self.headers['bonds'] += len(newbonds)
|
||
|
|
self.sections['Atoms'] = []
|
||
|
|
for atom in self.atoms + newatoms:
|
||
|
|
self.sections['Atoms'].append(atomline(atom))
|
||
|
|
for bond in newbonds:
|
||
|
|
self.sections['Bonds'].append(bondline(bond))
|
||
|
|
if 'Velocities' in self.sections:
|
||
|
|
self.sections['Velocities'] = []
|
||
|
|
for atom in self.atoms + newatoms:
|
||
|
|
self.sections['Velocities'].append(velline(atom))
|
||
|
|
|
||
|
|
# update list of atom IDs
|
||
|
|
for att in newattypes:
|
||
|
|
self.atomtypes.append(att)
|
||
|
|
|
||
|
|
|
||
|
|
def extract_pol(self, drude):
|
||
|
|
"""extract atom, drude, bonds info from polarizable data"""
|
||
|
|
|
||
|
|
# extract atom IDs
|
||
|
|
for line in self.sections['Masses']:
|
||
|
|
tok = line.split()
|
||
|
|
atomtype = {}
|
||
|
|
atomtype['id'] = int(tok[0])
|
||
|
|
atomtype['m'] = float(tok[1])
|
||
|
|
if len(tok) >= 4:
|
||
|
|
atomtype['type'] = tok[3]
|
||
|
|
atomtype['dflag'] = 'n'
|
||
|
|
if tok[-1] == "DC":
|
||
|
|
atomtype['dflag'] = 'c'
|
||
|
|
elif tok[-1] == "DP":
|
||
|
|
atomtype['dflag'] = 'd'
|
||
|
|
print(atomtype['dflag'])
|
||
|
|
else:
|
||
|
|
raise RuntimeError("comments in Masses section required "\
|
||
|
|
"to identify cores (DC) and Drudes (DP)")
|
||
|
|
self.atomtypes.append(atomtype)
|
||
|
|
|
||
|
|
# extract bond type data
|
||
|
|
for line in self.sections['Bond Coeffs']:
|
||
|
|
tok = line.split()
|
||
|
|
bondtype = {}
|
||
|
|
bondtype['id'] = int(tok[0])
|
||
|
|
bondtype['k'] = float(tok[1])
|
||
|
|
bondtype['r0'] = float(tok[2])
|
||
|
|
bondtype['note'] = ''.join([s + ' ' for s in tok[3:]]).strip()
|
||
|
|
self.bondtypes.append(bondtype)
|
||
|
|
|
||
|
|
# extract atom registers
|
||
|
|
for line in self.sections['Atoms']:
|
||
|
|
tok = line.split()
|
||
|
|
atom = {}
|
||
|
|
atom['n'] = int(tok[0])
|
||
|
|
atom['mol'] = int(tok[1])
|
||
|
|
atom['id'] = int(tok[2])
|
||
|
|
atom['q'] = float(tok[3])
|
||
|
|
atom['x'] = float(tok[4])
|
||
|
|
atom['y'] = float(tok[5])
|
||
|
|
atom['z'] = float(tok[6])
|
||
|
|
# atom['note'] = ''.join([s + ' ' for s in tok[7:-1]]).strip()
|
||
|
|
if tok[-1] == 'DC':
|
||
|
|
atom['note'] = ' '.join(tok[7:-1])
|
||
|
|
else:
|
||
|
|
atom['note'] = ' '.join(tok[7:])
|
||
|
|
self.atoms.append(atom)
|
||
|
|
self.idmap[atom['n']] = atom
|
||
|
|
|
||
|
|
if 'Velocities' in self.sections:
|
||
|
|
for line in self.sections['Velocities']:
|
||
|
|
tok = line.split()
|
||
|
|
atom = self.idmap[int(tok[0])]
|
||
|
|
atom['vx'] = float(tok[1])
|
||
|
|
atom['vy'] = float(tok[2])
|
||
|
|
atom['vz'] = float(tok[3])
|
||
|
|
|
||
|
|
# extract bond data
|
||
|
|
for line in self.sections['Bonds']:
|
||
|
|
tok = line.split()
|
||
|
|
bond = {}
|
||
|
|
bond['n'] = int(tok[0])
|
||
|
|
bond['id'] = int(tok[1])
|
||
|
|
bond['i'] = int(tok[2])
|
||
|
|
bond['j'] = int(tok[3])
|
||
|
|
bond['note'] = ''.join([s + ' ' for s in tok[4:]]).strip()
|
||
|
|
self.bonds.append(bond)
|
||
|
|
|
||
|
|
if 'Velocities' in self.sections:
|
||
|
|
for line in self.sections['Velocities']:
|
||
|
|
tok = line.split()
|
||
|
|
atom = self.idmap[int(tok[0])]
|
||
|
|
atom['vx'] = float(tok[1])
|
||
|
|
atom['vy'] = float(tok[2])
|
||
|
|
atom['vz'] = float(tok[3])
|
||
|
|
|
||
|
|
|
||
|
|
def depolarize(self, drude):
|
||
|
|
"""remove Drude particles"""
|
||
|
|
|
||
|
|
self.extract_pol(drude)
|
||
|
|
|
||
|
|
atom_tp_map = {}
|
||
|
|
bond_tp_map = {}
|
||
|
|
atom_map = {}
|
||
|
|
bond_map = {}
|
||
|
|
q = {}
|
||
|
|
atom_tp = {}
|
||
|
|
m = {}
|
||
|
|
|
||
|
|
for att in self.atomtypes:
|
||
|
|
if att['dflag'] != 'd':
|
||
|
|
atom_tp_map[att['id']] = len(atom_tp_map) + 1
|
||
|
|
m[att['id']] = att['m']
|
||
|
|
print(atom_tp_map)
|
||
|
|
for atom in self.atoms:
|
||
|
|
if atom['id'] in atom_tp_map:
|
||
|
|
atom_map[atom['n']] = len(atom_map) + 1
|
||
|
|
q[atom['n']] = atom['q']
|
||
|
|
atom_tp[atom['n']] = atom['id']
|
||
|
|
for bond in self.bonds:
|
||
|
|
if bond['i'] in atom_map and bond['j'] in atom_map:
|
||
|
|
bond_map[bond['n']] = len(bond_map) + 1
|
||
|
|
if bond['id'] not in bond_tp_map:
|
||
|
|
bond_tp_map[bond['id']] = len(bond_tp_map) + 1
|
||
|
|
else:
|
||
|
|
if bond['i'] in atom_map:
|
||
|
|
q[bond['i']] += q[bond['j']]
|
||
|
|
if atom_tp[bond['j']] in m:
|
||
|
|
m[atom_tp[bond['i']]] += m.pop(atom_tp[bond['j']])
|
||
|
|
else:
|
||
|
|
q[bond['j']] += q[bond['i']]
|
||
|
|
if atom_tp[bond['i']] in m:
|
||
|
|
m[atom_tp[bond['j']]] += m.pop(atom_tp[bond['i']])
|
||
|
|
|
||
|
|
size = len(self.atomtypes)
|
||
|
|
for iatom_tp in reversed(range(size)):
|
||
|
|
att = self.atomtypes[iatom_tp]
|
||
|
|
if att['id'] not in atom_tp_map:
|
||
|
|
del self.atomtypes[iatom_tp]
|
||
|
|
else:
|
||
|
|
att['m'] = m[att['id']]
|
||
|
|
att['id'] = atom_tp_map[att['id']]
|
||
|
|
|
||
|
|
size = len(self.bondtypes)
|
||
|
|
for ibond_tp in reversed(range(size)):
|
||
|
|
bdt = self.bondtypes[ibond_tp]
|
||
|
|
if bdt['id'] not in bond_tp_map:
|
||
|
|
del self.bondtypes[ibond_tp]
|
||
|
|
else:
|
||
|
|
bdt['id'] = bond_tp_map[bdt['id']]
|
||
|
|
|
||
|
|
size = len(self.atoms)
|
||
|
|
for iatom in reversed(range(size)):
|
||
|
|
atom = self.atoms[iatom]
|
||
|
|
if atom['n'] not in atom_map:
|
||
|
|
del self.atoms[iatom]
|
||
|
|
else:
|
||
|
|
atom['q'] = q[atom['n']]
|
||
|
|
atom['n'] = atom_map[atom['n']]
|
||
|
|
atom['id'] = atom_tp_map[atom['id']]
|
||
|
|
|
||
|
|
size = len(self.bonds)
|
||
|
|
for ibond in reversed(range(size)):
|
||
|
|
bond = self.bonds[ibond]
|
||
|
|
if bond['n'] not in bond_map:
|
||
|
|
del self.bonds[ibond]
|
||
|
|
else:
|
||
|
|
bond['n'] = bond_map[bond['n']]
|
||
|
|
bond['id'] = bond_tp_map[bond['id']]
|
||
|
|
bond['i'] = atom_map[bond['i']]
|
||
|
|
bond['j'] = atom_map[bond['j']]
|
||
|
|
|
||
|
|
self.sections['Atoms'] = []
|
||
|
|
for atom in self.atoms:
|
||
|
|
self.sections['Atoms'].append(atomline(atom))
|
||
|
|
self.headers['atoms'] = len(self.atoms)
|
||
|
|
self.sections['Masses'] = []
|
||
|
|
for att in self.atomtypes:
|
||
|
|
self.sections['Masses'].append(massline(att))
|
||
|
|
self.headers['atom types'] = len(self.atomtypes)
|
||
|
|
self.sections['Bonds'] = []
|
||
|
|
for bond in self.bonds:
|
||
|
|
self.sections['Bonds'].append(bondline(bond))
|
||
|
|
self.headers['bonds'] = len(self.bonds)
|
||
|
|
self.sections['Bond Coeffs'] = []
|
||
|
|
for bdt in self.bondtypes:
|
||
|
|
self.sections['Bond Coeffs'].append(bdtline(bdt))
|
||
|
|
self.headers['bond types'] = len(self.bondtypes)
|
||
|
|
|
||
|
|
if 'Velocities' in self.sections:
|
||
|
|
self.sections['Velocities'] = []
|
||
|
|
for atom in self.atoms:
|
||
|
|
self.sections['Velocities'].append(velline(atom))
|
||
|
|
|
||
|
|
def lmpscript(self, drude, outfile, thole = 2.6, cutoff = 12.0):
|
||
|
|
"""print lines for input script, including pair_style thole"""
|
||
|
|
|
||
|
|
dfound = False
|
||
|
|
for att in self.atomtypes:
|
||
|
|
if att['dflag'] == 'd':
|
||
|
|
dfound = True
|
||
|
|
break
|
||
|
|
if not dfound:
|
||
|
|
print("# No polarizable atoms found.")
|
||
|
|
return
|
||
|
|
|
||
|
|
print("# Commands to include in the LAMMPS input script\n")
|
||
|
|
|
||
|
|
print("# adapt the pair_style line as needed")
|
||
|
|
print("pair_style hybrid/overlay ... coul/long {0:.1f} "\
|
||
|
|
"thole {1:.3f} {0:.1f}\n".format(cutoff, thole))
|
||
|
|
|
||
|
|
print("read_data {0}\n".format(outfile))
|
||
|
|
|
||
|
|
print("# add interactions between any atoms and Drude particles")
|
||
|
|
print("pair_coeff * {0}* coul/long".format(att['id']))
|
||
|
|
|
||
|
|
# Thole parameters for I,J pairs
|
||
|
|
print("# add Thole screening if more than 1 Drude per molecule")
|
||
|
|
ifound = False
|
||
|
|
for atti in self.atomtypes:
|
||
|
|
itype = atti['type'].split()[0]
|
||
|
|
for ddt in drude.types:
|
||
|
|
dtype = ddt['type'].split()[0]
|
||
|
|
if dtype == itype:
|
||
|
|
alphai = ddt['alpha']
|
||
|
|
tholei = ddt['thole']
|
||
|
|
ifound = True
|
||
|
|
break
|
||
|
|
jfound = False
|
||
|
|
for attj in self.atomtypes:
|
||
|
|
if attj['id'] < atti['id']:
|
||
|
|
continue
|
||
|
|
jtype = attj['type'].split()[0]
|
||
|
|
for ddt in drude.types:
|
||
|
|
dtype = ddt['type'].split()[0]
|
||
|
|
if dtype == jtype:
|
||
|
|
alphaj = ddt['alpha']
|
||
|
|
tholej = ddt['thole']
|
||
|
|
jfound = True
|
||
|
|
break
|
||
|
|
if ifound and jfound:
|
||
|
|
alphaij = (alphai * alphaj)**0.5
|
||
|
|
tholeij = (tholei + tholej) / 2.0
|
||
|
|
print("pair_coeff {0:4} {1:4} thole {2:7.3f} "\
|
||
|
|
"{3:7.3f}".format(atti['id'], attj['id'],
|
||
|
|
alphaij, tholeij))
|
||
|
|
jfound = False
|
||
|
|
ifound = False
|
||
|
|
print("")
|
||
|
|
|
||
|
|
print("# atom groups convenient for thermostats (see package "
|
||
|
|
"documentation), etc.")
|
||
|
|
gatoms = gcores = gdrudes = ""
|
||
|
|
for att in self.atomtypes:
|
||
|
|
if att['dflag'] != 'd':
|
||
|
|
gatoms += " {0}".format(att['id'])
|
||
|
|
if att['dflag'] == 'c':
|
||
|
|
gcores += " {0}".format(att['id'])
|
||
|
|
if att['dflag'] == 'd':
|
||
|
|
gdrudes += " {0}".format(att['id'])
|
||
|
|
print("group ATOMS type" + gatoms)
|
||
|
|
print("group CORES type" + gcores)
|
||
|
|
print("group DRUDES type" + gdrudes)
|
||
|
|
print("")
|
||
|
|
|
||
|
|
print("# flag for each atom type: [C]ore, [D]rude, [N]on-polarizable")
|
||
|
|
drudetypes = ""
|
||
|
|
for att in self.atomtypes:
|
||
|
|
drudetypes += " {0}".format(att['dflag'].upper())
|
||
|
|
print("fix DRUDE all drude" + drudetypes)
|
||
|
|
print("")
|
||
|
|
|
||
|
|
print("# ATTENTION!")
|
||
|
|
print("# * special_bonds may need 'extra' keyword, LAMMPS will exit "
|
||
|
|
"with a message")
|
||
|
|
print("# * give all I<=J pair interactions, no mixing")
|
||
|
|
print("# * if using fix shake the group-ID must not include "
|
||
|
|
"Drude particles")
|
||
|
|
print("# use group ATOMS for example")
|
||
|
|
|
||
|
|
# --------------------------------------
|
||
|
|
|
||
|
|
kcal = 4.184 # kJ
|
||
|
|
eV = 96.485 # kJ/mol
|
||
|
|
fpe0 = 0.000719756 # (4 Pi eps0) in e^2/(kJ/mol A)
|
||
|
|
|
||
|
|
|
||
|
|
class Drude(object):
|
||
|
|
"""specification of drude oscillator types"""
|
||
|
|
|
||
|
|
def __init__(self, drudefile, polar = '', positive = False, metal = False):
|
||
|
|
self.types = []
|
||
|
|
with open(drudefile, "r") as f:
|
||
|
|
for line in f:
|
||
|
|
line = line.strip()
|
||
|
|
if line.startswith('#') or len(line) == 0:
|
||
|
|
continue
|
||
|
|
tok = line.split()
|
||
|
|
drude = {}
|
||
|
|
drude['type'] = tok[0]
|
||
|
|
drude['dm'] = float(tok[1])
|
||
|
|
dq = float(tok[2])
|
||
|
|
k = float(tok[3])
|
||
|
|
drude['alpha'] = alpha = float(tok[4])
|
||
|
|
drude['thole'] = float(tok[5])
|
||
|
|
|
||
|
|
if polar == 'q':
|
||
|
|
dq = (fpe0 * k * alpha)**0.5
|
||
|
|
elif polar == 'k':
|
||
|
|
k = dq*dq / (fpe0 * alpha)
|
||
|
|
|
||
|
|
if positive:
|
||
|
|
drude['dq'] = abs(dq)
|
||
|
|
else:
|
||
|
|
drude['dq'] = -abs(dq)
|
||
|
|
|
||
|
|
if metal:
|
||
|
|
drude['k'] = k / (2.0 * eV)
|
||
|
|
else:
|
||
|
|
drude['k'] = k / (2.0 * kcal)
|
||
|
|
|
||
|
|
self.types.append(drude)
|
||
|
|
|
||
|
|
|
||
|
|
# --------------------------------------
|
||
|
|
|
||
|
|
def main():
|
||
|
|
parser = argparse.ArgumentParser(description = usage,
|
||
|
|
formatter_class = argparse.RawTextHelpFormatter)
|
||
|
|
parser.add_argument('-f', '--ffdrude', default = 'drude.dff',
|
||
|
|
help = 'Drude parameter file (default: drude.dff)')
|
||
|
|
parser.add_argument('-t', '--thole', type = float, default = 2.6,
|
||
|
|
help = 'Thole damping parameter (default: 2.6)')
|
||
|
|
parser.add_argument('-c', '--cutoff', type = float, default = 12.0,
|
||
|
|
help = 'distance cutoff/A (default: 12.0)')
|
||
|
|
parser.add_argument('-q', '--qcalc', action = 'store_true',
|
||
|
|
help = 'Drude charges calculated from polarisability '\
|
||
|
|
'(default: q value from parameter file)')
|
||
|
|
parser.add_argument('-k', '--kcalc', action = 'store_true',
|
||
|
|
help = 'Drude force constants calculated from '\
|
||
|
|
'polarisability (default: k value from parameter file)')
|
||
|
|
parser.add_argument('-p', '--positive', action = 'store_true',
|
||
|
|
help = 'Drude particles have positive charge '\
|
||
|
|
'(default: negative charge)')
|
||
|
|
parser.add_argument('-m', '--metal', action = 'store_true',
|
||
|
|
help = 'LAMMPS metal units (default: real units)')
|
||
|
|
parser.add_argument('-d', '--depolarize', action = 'store_true',
|
||
|
|
help = 'remove Drude dipole polarization from '\
|
||
|
|
'LAMMPS data file')
|
||
|
|
parser.add_argument('infile', help = 'input LAMMPS data file')
|
||
|
|
parser.add_argument('outfile', help = 'output LAMMPS data file')
|
||
|
|
args = parser.parse_args()
|
||
|
|
|
||
|
|
if args.qcalc:
|
||
|
|
polar = 'q'
|
||
|
|
elif args.kcalc:
|
||
|
|
polar = 'p'
|
||
|
|
else:
|
||
|
|
polar = ''
|
||
|
|
|
||
|
|
data = Data(args.infile)
|
||
|
|
drude = Drude(args.ffdrude, polar, args.positive, args.metal)
|
||
|
|
if not args.depolarize:
|
||
|
|
data.polarize(drude)
|
||
|
|
data.lmpscript(drude, args.outfile, args.thole, args.cutoff)
|
||
|
|
else:
|
||
|
|
data.depolarize(drude)
|
||
|
|
data.write(args.outfile)
|
||
|
|
|
||
|
|
if __name__ == '__main__':
|
||
|
|
main()
|