mirror of
https://github.com/ParticulateFlow/LPP.git
synced 2025-12-08 06:37:46 +00:00
269 lines
8.3 KiB
Python
269 lines
8.3 KiB
Python
# Pizza.py toolkit, www.cs.sandia.gov/~sjplimp/pizza.html
|
|
# Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
|
|
#
|
|
# Copyright (2005) Sandia Corporation. Under the terms of Contract
|
|
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
|
# certain rights in this software. This software is distributed under
|
|
# the GNU General Public License.
|
|
|
|
# pair tool
|
|
|
|
# Imports and external programs
|
|
|
|
from __future__ import print_function, absolute_import
|
|
from math import sqrt
|
|
|
|
oneline = "Compute LAMMPS pairwise energies"
|
|
|
|
docstr = """
|
|
p = pair("lj/charmm/coul/charmm") create pair object for specific pair style
|
|
|
|
available styles: lj/cut, lj/cut/coul/cut, lj/charmm/coul/charmm
|
|
|
|
p.coeff(d) extract pairwise coeffs from data object
|
|
p.init(cut1,cut2,...) setup based on coeffs and cutoffs
|
|
|
|
init args are specific to pair style:
|
|
lj/cut = cutlj
|
|
lj/cut/coul/cut = cutlj,cut_coul (cut_coul optional)
|
|
lj/charmm/coul/charmm = cutlj_inner,cutlj,cutcoul_inner,cut_coul
|
|
(last 2 optional)
|
|
|
|
e_vdwl,e_coul = p.single(rsq,itype,jtype,q1,q2,...) compute LJ/Coul energy
|
|
|
|
pairwise energy between 2 atoms at distance rsq with their attributes
|
|
args are specific to pair style:
|
|
lj/cut = rsq,itype,jtype
|
|
lj/cut/coul/cut = rsq,itype,jtype,q1,q2
|
|
lj/charmm/coul/charmm = rsq,itype,jtype,q1,q2
|
|
"""
|
|
|
|
# History
|
|
# 8/05, Steve Plimpton and Paul Crozier (SNL): original version
|
|
# 9/05, Paul Crozier (SNL): added lj/cut and lj/cut/coul/cut
|
|
|
|
# ToDo list
|
|
|
|
# Variables
|
|
|
|
# Class definition
|
|
|
|
class pair:
|
|
|
|
# --------------------------------------------------------------------
|
|
|
|
def __init__(self,style):
|
|
|
|
# map generic methods to style-specific methods
|
|
|
|
if style == "lj/cut":
|
|
self.coeff_func = self.coeff_lj_cut
|
|
self.init_func = self.init_lj_cut
|
|
self.single_func = self.single_lj_cut
|
|
elif style == "lj/cut/coul/cut":
|
|
self.coeff_func = self.coeff_lj_cut_coul_cut
|
|
self.init_func = self.init_lj_cut_coul_cut
|
|
self.single_func = self.single_lj_cut_coul_cut
|
|
elif style == "lj/charmm/coul/charmm":
|
|
self.coeff_func = self.coeff_lj_charmm_coul_charmm
|
|
self.init_func = self.init_lj_charmm_coul_charmm
|
|
self.single_func = self.single_lj_charmm_coul_charmm
|
|
else:
|
|
raise Exception("this pair style not yet supported")
|
|
|
|
# --------------------------------------------------------------------
|
|
# generic coeff method
|
|
|
|
def coeff(self,data):
|
|
self.coeff_func(data)
|
|
|
|
# --------------------------------------------------------------------
|
|
# generic init method, as many args as needed
|
|
|
|
def init(self,*list):
|
|
self.init_func(list)
|
|
|
|
# --------------------------------------------------------------------
|
|
# generic single method, as many args as needed
|
|
|
|
def single(self,*list):
|
|
return self.single_func(list)
|
|
|
|
# --------------------------------------------------------------------
|
|
# --------------------------------------------------------------------
|
|
# lj/cut methods
|
|
|
|
def coeff_lj_cut(self,data):
|
|
epsilon = data.get("Pair Coeffs",2)
|
|
sigma = data.get("Pair Coeffs",3)
|
|
ntypes = len(epsilon)
|
|
|
|
self.lj3 = []
|
|
self.lj4 = []
|
|
for i in range(ntypes):
|
|
self.lj3.append(ntypes * [0])
|
|
self.lj4.append(ntypes * [0])
|
|
for j in range(ntypes):
|
|
epsilon_ij = sqrt(epsilon[i]*epsilon[j])
|
|
sigma_ij = sqrt(sigma[i]*sigma[j])
|
|
self.lj3[i][j] = 4.0 * epsilon_ij * pow(sigma_ij,12.0);
|
|
self.lj4[i][j] = 4.0 * epsilon_ij * pow(sigma_ij,6.0);
|
|
|
|
# --------------------------------------------------------------------
|
|
# args = cutlj
|
|
|
|
def init_lj_cut(self,list):
|
|
cut_lj = list[0]
|
|
self.cut_ljsq = cut_lj*cut_lj
|
|
|
|
# --------------------------------------------------------------------
|
|
# args = rsq,itype,jtype
|
|
|
|
def single_lj_cut(self,list):
|
|
rsq = list[0]
|
|
itype = list[1]
|
|
jtype = list[2]
|
|
|
|
r2inv = 1.0/rsq
|
|
|
|
if rsq < self.cut_ljsq:
|
|
r6inv = r2inv*r2inv*r2inv
|
|
eng_vdwl = r6inv*(self.lj3[itype][jtype]*r6inv-self.lj4[itype][jtype])
|
|
else: eng_vdwl = 0.0
|
|
|
|
return eng_vdwl
|
|
|
|
# --------------------------------------------------------------------
|
|
# --------------------------------------------------------------------
|
|
# lj/cut/coul/cut methods
|
|
|
|
def coeff_lj_cut_coul_cut(self,data):
|
|
epsilon = data.get("Pair Coeffs",2)
|
|
sigma = data.get("Pair Coeffs",3)
|
|
ntypes = len(epsilon)
|
|
|
|
self.lj3 = []
|
|
self.lj4 = []
|
|
for i in range(ntypes):
|
|
self.lj3.append(ntypes * [0])
|
|
self.lj4.append(ntypes * [0])
|
|
for j in range(ntypes):
|
|
epsilon_ij = sqrt(epsilon[i]*epsilon[j])
|
|
sigma_ij = sqrt(sigma[i]*sigma[j])
|
|
self.lj3[i][j] = 4.0 * epsilon_ij * pow(sigma_ij,12.0);
|
|
self.lj4[i][j] = 4.0 * epsilon_ij * pow(sigma_ij,6.0);
|
|
|
|
# --------------------------------------------------------------------
|
|
# args = cutlj, cut_coul (cut_coul optional)
|
|
|
|
def init_lj_cut_coul_cut(self,list):
|
|
self.qqr2e = 332.0636 # convert energy to kcal/mol
|
|
cut_lj = list[0]
|
|
self.cut_ljsq = cut_lj*cut_lj
|
|
|
|
if len(list) == 1: cut_coul = cut_lj
|
|
else: cut_coul = list[1]
|
|
self.cut_coulsq = cut_coul*cut_coul
|
|
|
|
# --------------------------------------------------------------------
|
|
# args = rsq,itype,jtype,q1,q2
|
|
|
|
def single_lj_cut_coul_cut(self,list):
|
|
rsq = list[0]
|
|
itype = list[1]
|
|
jtype = list[2]
|
|
q1 = list[3]
|
|
q2 = list[4]
|
|
|
|
r2inv = 1.0/rsq
|
|
|
|
if rsq < self.cut_coulsq: eng_coul = self.qqr2e * q1*q2*sqrt(r2inv)
|
|
else: eng_coul = 0.0
|
|
|
|
if rsq < self.cut_ljsq:
|
|
r6inv = r2inv*r2inv*r2inv
|
|
eng_vdwl = r6inv*(self.lj3[itype][jtype]*r6inv-self.lj4[itype][jtype])
|
|
else: eng_vdwl = 0.0
|
|
|
|
return eng_coul,eng_vdwl
|
|
|
|
# --------------------------------------------------------------------
|
|
# --------------------------------------------------------------------
|
|
# lj/charmm/coul/charmm methods
|
|
|
|
def coeff_lj_charmm_coul_charmm(self,data):
|
|
epsilon = data.get("Pair Coeffs",2)
|
|
sigma = data.get("Pair Coeffs",3)
|
|
ntypes = len(epsilon)
|
|
|
|
self.lj3 = []
|
|
self.lj4 = []
|
|
for i in range(ntypes):
|
|
self.lj3.append(ntypes * [0])
|
|
self.lj4.append(ntypes * [0])
|
|
for j in range(ntypes):
|
|
epsilon_ij = sqrt(epsilon[i]*epsilon[j])
|
|
sigma_ij = 0.5 * (sigma[i] + sigma[j])
|
|
self.lj3[i][j] = 4.0 * epsilon_ij * pow(sigma_ij,12.0);
|
|
self.lj4[i][j] = 4.0 * epsilon_ij * pow(sigma_ij,6.0);
|
|
|
|
# --------------------------------------------------------------------
|
|
# args = cutlj_inner,cutlj,cutcoul_inner,cut_coul (last 2 optional)
|
|
|
|
def init_lj_charmm_coul_charmm(self,list):
|
|
self.qqr2e = 332.0636 # convert energy to kcal/mol
|
|
cut_lj_inner = list[0]
|
|
cut_lj = list[1]
|
|
|
|
self.cut_lj_innersq = cut_lj_inner*cut_lj_inner
|
|
self.cut_ljsq = cut_lj*cut_lj
|
|
|
|
if len(list) == 2:
|
|
cut_coul_inner = cut_lj_inner
|
|
cut_coul = cut_lj
|
|
else:
|
|
cut_coul_inner = list[2]
|
|
cut_coul = list[3]
|
|
self.cut_coul_innersq = cut_coul_inner*cut_coul_inner
|
|
self.cut_coulsq = cut_coul*cut_coul
|
|
|
|
self.denom_lj = (self.cut_ljsq-self.cut_lj_innersq) * \
|
|
(self.cut_ljsq-self.cut_lj_innersq) * \
|
|
(self.cut_ljsq-self.cut_lj_innersq);
|
|
self.denom_coul = (self.cut_coulsq-self.cut_coul_innersq) * \
|
|
(self.cut_coulsq-self.cut_coul_innersq) * \
|
|
(self.cut_coulsq-self.cut_coul_innersq);
|
|
|
|
# --------------------------------------------------------------------
|
|
# args = rsq,itype,jtype,q1,q2
|
|
|
|
def single_lj_charmm_coul_charmm(self,list):
|
|
rsq = list[0]
|
|
itype = list[1]
|
|
jtype = list[2]
|
|
q1 = list[3]
|
|
q2 = list[4]
|
|
|
|
r2inv = 1.0/rsq
|
|
|
|
if rsq < self.cut_coulsq:
|
|
eng_coul = self.qqr2e * q1*q2*sqrt(r2inv)
|
|
if rsq > self.cut_coul_innersq:
|
|
switch1 = (self.cut_coulsq-rsq) * (self.cut_coulsq-rsq) * \
|
|
(self.cut_coulsq + 2.0*rsq - 3.0*self.cut_coul_innersq) / \
|
|
self.denom_coul
|
|
eng_coul *= switch1
|
|
else: eng_coul = 0.0
|
|
|
|
if rsq < self.cut_ljsq:
|
|
r6inv = r2inv*r2inv*r2inv
|
|
eng_vdwl = r6inv*(self.lj3[itype][jtype]*r6inv-self.lj4[itype][jtype])
|
|
if rsq > self.cut_lj_innersq:
|
|
switch1 = (self.cut_ljsq-rsq) * (self.cut_ljsq-rsq) * \
|
|
(self.cut_ljsq + 2.0*rsq - 3.0*self.cut_lj_innersq) / \
|
|
self.denom_lj
|
|
eng_vdwl *= switch1
|
|
else: eng_vdwl = 0.0
|
|
|
|
return eng_coul,eng_vdwl
|