mirror of https://github.com/lammps/lammps.git
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13730 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
3dfb446590
commit
d58afbcc48
|
|
@ -19,6 +19,7 @@ chain create a data file of bead-spring chains
|
|||
colvars post-process output of the fix colvars command
|
||||
createatoms generate lattices of atoms within a geometry
|
||||
data2xmovie convert a data file to a snapshot that xmovie can viz
|
||||
drude create Drude core/electron atom pairs in a data file
|
||||
eam_database one tool to generate EAM alloy potential files
|
||||
eam_generate 2nd tool to generate EAM alloy potential files
|
||||
eff scripts for working with the eFF (electron force field)
|
||||
|
|
|
|||
|
|
@ -0,0 +1,29 @@
|
|||
# drude.ff, version 2015/05/07
|
||||
# units: kJ/mol, A, deg
|
||||
# kforce is in the form k/2 r_D^2
|
||||
|
||||
# type m_D/u q_D/e k_D alpha/A3 thole
|
||||
# PSPC water, JCP 119(2003)3025
|
||||
Ow 0.8 2.0824 4184.0 1.440 2.6
|
||||
# SWM4-NDP water, CPL 418(2006)245
|
||||
ODw 0.8 -1.7164 4184.0 0.978 2.6
|
||||
# alkanes (CHARMM)
|
||||
C3H 0.4 -1.0 4184.0 2.051 2.6
|
||||
C2H 0.4 -1.0 4184.0 1.660 2.6
|
||||
C1H 0.4 -1.0 4184.0 1.302 2.6
|
||||
C4C 0.4 -1.0 4184.0 1.302 2.6
|
||||
# alcohols (CHARMM)
|
||||
#OH 0.4 -1.0 4184.0 1.028 2.6
|
||||
#CTO 0.4 -1.0 4184.0 1.000 2.6
|
||||
#C2O 0.4 -1.0 4184.0 1.400 2.6
|
||||
# ethanol Noskov JPCB 109(2005)6705
|
||||
OH 0.4 -1.0 4184.0 0.97772 2.6
|
||||
CTO 0.4 -1.0 4184.0 1.21757 2.6
|
||||
C2O 0.4 -1.0 4184.0 1.23758 2.6
|
||||
# aromatic (CHARMM)
|
||||
CA 0.4 -1.0 4184.0 1.62 2.6
|
||||
CG 0.4 -1.0 4184.0 1.62 2.6
|
||||
# general elements, Schwerdtfeger Table 2014
|
||||
S 0.8 -1.0 4184.0 2.90 2.6
|
||||
P 0.8 -1.0 4184.0 3.69 2.6
|
||||
Mo 0.8 -1.0 4184.0 10.74 2.6
|
||||
|
|
@ -0,0 +1,726 @@
|
|||
#!/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()
|
||||
|
|
@ -0,0 +1,70 @@
|
|||
# general elements, Schwerdtfeger Table 2014
|
||||
# units: kJ/mol, A, deg
|
||||
# kforce is in the form k/2 r_D^2
|
||||
|
||||
# type m_D/u q_D/e k_D alpha/A3 thole
|
||||
H 0.8 -1.0 4184.0 0.6679 2.6
|
||||
He 0.8 -1.0 4184.0 0.2051 2.6
|
||||
Li 0.8 -1.0 4184.0 24.3023 2.6
|
||||
Be 0.8 -1.0 4184.0 5.5880 2.6
|
||||
B 0.8 -1.0 4184.0 3.0422 2.6
|
||||
C 0.8 -1.0 4184.0 1.6686 2.6
|
||||
N 0.8 -1.0 4184.0 1.1262 2.6
|
||||
O 0.8 -1.0 4184.0 0.7765 2.6
|
||||
F 0.8 -1.0 4184.0 0.5483 2.6
|
||||
Ne 0.8 -1.0 4184.0 0.3957 2.6
|
||||
Na 0.8 -1.0 4184.0 24.1097 2.6
|
||||
Mg 0.8 -1.0 4184.0 10.4856 2.6
|
||||
Al 0.8 -1.0 4184.0 6.8165 2.6
|
||||
Si 0.8 -1.0 4184.0 5.5288 2.6
|
||||
P 0.8 -1.0 4184.0 3.6942 2.6
|
||||
S 0.8 -1.0 4184.0 2.8703 2.6
|
||||
Cl 0.8 -1.0 4184.0 2.1591 2.6
|
||||
Ar 0.8 -1.0 4184.0 1.6404 2.6
|
||||
K 0.8 -1.0 4184.0 43.0625 2.6
|
||||
Ca 0.8 -1.0 4184.0 25.0432 2.6
|
||||
Sc 0.8 -1.0 4184.0 21.0837 2.6
|
||||
Ti 0.8 -1.0 4184.0 16.9434 2.6
|
||||
V 0.8 -1.0 4184.0 14.4243 2.6
|
||||
Cr 0.8 -1.0 4184.0 11.6177 2.6
|
||||
Mn 0.8 -1.0 4184.0 9.8987 2.6
|
||||
Fe 0.8 -1.0 4184.0 9.2838 2.6
|
||||
Co 0.8 -1.0 4184.0 8.5517 2.6
|
||||
Ni 0.8 -1.0 4184.0 7.5722 2.6
|
||||
Cu 0.8 -1.0 4184.0 6.0311 2.6
|
||||
Zn 0.8 -1.0 4184.0 5.7496 2.6
|
||||
Ga 0.8 -1.0 4184.0 7.6167 2.6
|
||||
Ge 0.8 -1.0 4184.0 5.8429 2.6
|
||||
As 0.8 -1.0 4184.0 4.4159 2.6
|
||||
Se 0.8 -1.0 4184.0 3.8884 2.6
|
||||
Br 0.8 -1.0 4184.0 3.1163 2.6
|
||||
Kr 0.8 -1.0 4184.0 2.5303 2.6
|
||||
Rb 0.8 -1.0 4184.0 47.2413 2.6
|
||||
Sr 0.8 -1.0 4184.0 27.5624 2.6
|
||||
Y 0.8 -1.0 4184.0 22.6723 2.6
|
||||
Zr 0.8 -1.0 4184.0 17.9304 2.6
|
||||
Nb 0.8 -1.0 4184.0 15.7076 2.6
|
||||
Mo 0.8 -1.0 4184.0 10.7434 2.6
|
||||
Tc 0.8 -1.0 4184.0 11.9141 2.6
|
||||
Ru 0.8 -1.0 4184.0 9.6320 2.6
|
||||
Rh 0.8 -1.0 4184.0 8.5947 2.6
|
||||
Pd 0.8 -1.0 4184.0 4.7419 2.6
|
||||
Ag 0.8 -1.0 4184.0 5.4384 2.6
|
||||
Cd 0.8 -1.0 4184.0 7.3574 2.6
|
||||
In 0.8 -1.0 4184.0 10.1803 2.6
|
||||
Sn 0.8 -1.0 4184.0 6.2830 2.6
|
||||
Sb 0.8 -1.0 4184.0 6.3053 2.6
|
||||
Te 0.8 -1.0 4184.0 5.4828 2.6
|
||||
I 0.8 -1.0 4184.0 5.1272 2.6
|
||||
Xe 0.8 -1.0 4184.0 4.1218 2.6
|
||||
Cs 0.8 -1.0 4184.0 59.4221 2.6
|
||||
Ba 0.8 -1.0 4184.0 39.7135 2.6
|
||||
La 0.8 -1.0 4184.0 31.1188 2.6
|
||||
Ce 0.8 -1.0 4184.0 29.6369 2.6
|
||||
Yb 0.8 -1.0 4184.0 21.1311 2.6
|
||||
Pt 0.8 -1.0 4184.0 6.5201 2.6
|
||||
Au 0.8 -1.0 4184.0 4.1344 2.6
|
||||
Hg 0.8 -1.0 4184.0 5.0249 2.6
|
||||
Pb 0.8 -1.0 4184.0 6.9795 2.6
|
||||
Bi 0.8 -1.0 4184.0 7.8316 2.6
|
||||
U 0.8 -1.0 4184.0 20.3013 2.6
|
||||
Loading…
Reference in New Issue