#!/usr/bin/env python3
'''
.. code-block::
!---------------------------------------------------------------------------!
! initial_conditions: Module for generating initial conditions !
! Implementations by: Yi-Fan Hou & Pavlo O. Dral !
!---------------------------------------------------------------------------!
'''
import numpy as np
import warnings
from . import data
from . import stopper
from . import constants
from . import thermo
[docs]
def generate_initial_conditions(molecule=None, generation_method=None, number_of_initial_conditions=1,
file_with_initial_xyz_coordinates=None, file_with_initial_xyz_velocities=None,
eliminate_angular_momentum=True,
degrees_of_freedom=None,
initial_temperature=None, initial_kinetic_energy=None):
'''
Generate initial conditions
Arguments:
molecule (:class:`data.molecule`): Molecule with necessary information
generation_method (str): Initial condition generation method
number_of_initial_conditions (int): Number of initial conditions to generate, 1 by default
file_with_initial_xyz_coordinates (str): File with initial xyz coordinates, only valid for ``generation_method='user-defined'``
file_with_initial_xyz_velocities (str): File with initial xyz velocities, only valid for ``generation_method='user-defined'``
eliminate_angular_momentum (bool): Remove angular momentum from velocities, valid for ``generation_method='random'``
degrees_of_freedom (int): Degrees of freedom of the molecule, by default remove translational and rotational degrees of freedom. It can be a negative value, which means that some value is subtracted from 3*Natoms
initial_temperature (float): Initial temperature in Kelvin, control random initial velocities
initial_kinetic_energy (float): Initial temperature in Hartree, control random initial velocities
.. table::
:align: center
============================= =============================================
generation_method description
============================= =============================================
``'user-defined'`` (default) Use user-defined initial conditions
``'random'`` Generate random velocities
============================= =============================================
.. ``'wigner'`` Use Wigner sampling
.. ``'maxwell-boltzmann'`` Randomly generate initial velocities from Maxwell-Boltzmann distribution
Returns:
A molecular database (:class:`ml.data.molecular_database`) with ``number_of_initial_conditions`` initial conditions
Examples:
.. code-block:: python
# Use user-defined initial conditions
init_cond_db = ml.generate_initial_conditions(molecule = mol,
generation_method = 'user-defined',
file_with_initial_xyz_coordinates = 'ethanol.xyz',
file_with_initial_xyz_velocities = 'ethanol.vxyz',
number_of_initial_conditions = 1)
# Generate random velocities
init_cond_db = ml.generate_initial_conditions(molecule = mol,
generation_method = 'random',
initial_temperature = 300,
number_of_initial_conditions = 1)
.. # Use Wigner sampling
.. init_cond_db = ml.generate_initial_conditions(molecule = mol,
.. generation_method = 'wigner',
.. number_of_initial_conditions = 1)
.. .. note::
.. Wigner sampling needs Hessian matrix. You can use ``ml.models.methods.predict(molecule=mol,calculate_hessian=True)`` to get Hessian matrix.
'''
if initial_temperature != None and initial_kinetic_energy != None:
stopper.stopMLatom('Cannot use initial_temperature and initial_kinetic_energy at the same time')
if initial_temperature == None and initial_kinetic_energy == None:
initial_temperature = 300
if degrees_of_freedom != None:
if degrees_of_freedom <= 0:
degrees_of_freedom = 3 * len(molecule.atoms) + degrees_of_freedom
else:
if eliminate_angular_momentum:
degrees_of_freedom = max(1, 3 * len(molecule.atoms) - 6)
else:
degrees_of_freedom = max(1, 3 * len(molecule.atoms) - 3)
init_cond_db = data.molecular_database()
Natoms = len(molecule.atoms)
if generation_method == 'user-defined':
init_cond_db.read_from_xyz_file(file_with_initial_xyz_coordinates)
init_cond_db.add_xyz_vectorial_properties_from_file(file_with_initial_xyz_velocities, xyz_vectorial_property='xyz_velocities')
elif generation_method == 'random':
for irepeat in range(number_of_initial_conditions):
new_molecule = molecule.copy(atomic_labels = ['xyz_coordinates'], molecular_labels = [])
velocities = generate_random_velocities(new_molecule,eliminate_angular_momentum,degrees_of_freedom,temp=initial_temperature,ekin=initial_kinetic_energy)
for iatom in range(Natoms):
new_molecule.atoms[iatom].xyz_velocities = velocities[iatom]
init_cond_db.molecules.append(new_molecule)
elif generation_method.casefold() == 'maxwell-boltzmann':
init_cond_db.read_from_numpy(coordinates=np.repeat([molecule.xyz_coordinates], number_of_initial_conditions, axis=0), species=np.repeat([molecule.element_symbols], number_of_initial_conditions, axis=0))
velocities = np.random.randn(*init_cond_db.xyz_coordinates.shape) * np.sqrt(initial_temperature * constants.kB / (init_cond_db.nuclear_masses / 1000 / constants.Avogadro_constant))[...,np.newaxis] / 1E5
init_cond_db.add_xyz_vectorial_properties(velocities, xyz_vectorial_property='xyz_velocities')
elif generation_method == 'wigner':
coordinates_all, velocities_all, linear = wignersample(number_of_initial_conditions,molecule)
mass_ = np.array([each.nuclear_mass for each in molecule.atoms])
mass = mass_.reshape(Natoms,1)
total_mass = np.sum(mass_)
for irepeat in range(number_of_initial_conditions):
if eliminate_angular_momentum:
velocities_all[irepeat] = getridofang(coordinates_all[irepeat],velocities_all[irepeat],mass_)
v_cm = sum(velocities_all[irepeat]*mass)/total_mass
velocities_all[irepeat] -= v_cm
for irepeat in range(number_of_initial_conditions):
new_molecule = molecule.copy(atomic_labels = [], molecular_labels = [])
for iatom in range(Natoms):
new_molecule.atoms[iatom].xyz_coordinates = coordinates_all[irepeat][iatom]
new_molecule.atoms[iatom].xyz_velocities = velocities_all[irepeat][iatom]
init_cond_db.molecules.append(new_molecule)
#print("wigner")
#init_cond_db.molecules.append(molecule) # Replace with actual new initial conditions
return init_cond_db
def read_velocities_from_file(filename):
velocities = []
with open(filename, 'r') as fxyz:
nlines = 0
natoms = 0
for line in fxyz:
nlines += 1
if nlines == 1:
natoms = int(line)
elif nlines > 2 and nlines <= 2 + natoms:
yy = line.split()
velocities.append([float(xx) for xx in yy[-3:]])
if nlines == 2 + natoms:
break
return velocities
def generate_random_velocities(molecule,noang,dof,temp=None,ekin=None):
np.random.seed()
Natoms = len(molecule.atoms)
randnum = np.random.randn(Natoms,3)
coord = np.array([each.xyz_coordinates for each in molecule.atoms])
mass_ = np.array([each.nuclear_mass for each in molecule.atoms])
mass = np.array(mass_).reshape(Natoms,1)
if temp != None:
# kb = 1.380649E-23 # Unit: J/K
# hartree = 27.21070 * 1.602176565E-19 # Unit: J/Hartree
# kb = kb / hartree # Unit: Hartree/K
kb = constants.kB_in_Hartree # Unit: Hartree/K
kinetic_energy = dof/2.0 * kb * temp
else:
kinetic_energy = ekin
linearity = molecule.is_it_linear()
# Eliminate total angular momentum
if noang:
if linearity:
rand_velocity = generate_random_velocities_for_linear_molecule(molecule)
else:
#rand_velocity = randnum / np.sqrt(mass*1822)
rand_velocity = randnum / np.sqrt(mass*constants.ram2au)
rand_velocity = getridofang(coord,rand_velocity,mass_)
else:
rand_velocity = randnum / np.sqrt(mass*constants.ram2au)
#print(rand_velocity)
# Raise warning if degrees of freedom are not compatible to the linear molecule
if linearity:
if dof != 3*Natoms-5:
#warnings.warn_explicit(f'Linear molecule detected, but degrees of freedom used is {dof} instead of {3*Natoms-5}')
print(f'WARNING: Linear molecule detected, but degrees of freedom used is {dof} instead of {3*Natoms-5}')
# Eliminate total linear momentum
total_mass = sum(mass)[0]
v_cm = sum(rand_velocity*mass)/total_mass
rand_velocity = rand_velocity - v_cm
#rand_energy = np.sum((rand_velocity**2) * (mass*1822)) / 2.0
rand_energy = np.sum((rand_velocity**2) * (mass*constants.ram2au)) / 2.0
ratio = rand_energy / kinetic_energy
velocity = rand_velocity / np.sqrt(ratio) # Unit: a.u.
# bohr = 0.5291772083 # 1 Bohr = 0.5291772083 Angstrom
# time = 2.4189E-2 # Time Unit: 1 a.u. = 2.4189E-17 s = 2.4189E-2 fs
# velocity = velocity * bohr / time # Unit: Angstrom / fs
velocity = velocity * constants.Bohr2Angstrom / constants.au2fs # Unit: Angstrom / fs
return velocity
def getridofang(coord,vel,mass):
omega = getAngV(coord,vel,mass)
coord_ = coord - getCoM(coord,mass)
vel = vel - np.cross(omega,coord)
return vel
def getAngV(xyz,v,m,center=None):
L=getAngM(xyz,v,m,center)
I=getMomentOfInertiaTensor(xyz,m,center)
omega=np.linalg.inv(I).dot(L)
return omega
def getCoM(xyz,m=None):
if m is None:
m=np.ones(xyz.shape[-2])
return np.sum(xyz*m[:,np.newaxis],axis=-2)/np.sum(m)
def getAngM(xyz,v,m,center=None):
if center is None:
centered=xyz-getCoM(xyz,m)
else:
centered=xyz-center
L=np.sum(m[:,np.newaxis]*np.cross(centered,v),axis=0)
return L
def getMomentOfInertiaTensor(xyz,m,center=None):
if center is None:
center=getCoM(xyz,m)
centered=xyz-center
I=np.zeros((3,3))
for i in range(3):
for j in range(3):
for k in range(len(m)):
I[i,j]+=m[k]*(np.sum(centered[k]**2)-centered[k,i]*centered[k,j]) if i==j else m[k]*(-centered[k,i]*centered[k,j])
return I
def readGaussianNM(filename):
freqs = []
ele = []
with open(filename,'r') as f:
lines = f.readlines()
freqind = []
for iline in range(len(lines)):
if 'Frequencies' in lines[iline]:
freqind.append(iline)
# Get the number of normal modes
Nnm = 0
for eachind in freqind:
Nnm += len(lines[eachind].strip().split()) - 2
# Get the number of atoms
Natoms = 0
indtemp = freqind[0]
ii = 1
while (True):
linesplit = lines[indtemp+ii]#.strip().split()
ii+=1
if 'Atom' in linesplit:
break
Ncols = 3 * (len(lines[indtemp].strip().split()) - 2) + 2
indtemp += ii
ii=0
while True:
if len(lines[indtemp+ii].strip().split()) == Ncols:
ii+=1
else:
break
Natoms = ii
# Read atom type
for iatom in range(Natoms):
ele.append(eval(lines[indtemp+iatom].strip().split()[1]))
# Read normal modes
nmflag = 0
ntriv = 3*Natoms - Nnm
if ntriv == 6:
linear = 0
elif ntriv == 5:
linear = 1
#print(linear)
normalmodes = np.zeros((3*Natoms,Natoms,3))
for itriv in range(ntriv):
normalmodes[itriv] = 1.0 / np.sqrt(3*Natoms)
for eachind_ in freqind:
#print(lines[eachind_])
Ncolnm = len(lines[eachind_].strip().split()) - 2
for icolnm in range(Ncolnm):
freqs.append(eval(lines[eachind_].strip().split()[2+icolnm]))
#print(freqs)
eachind = eachind_ + 5
for iatom in range(Natoms):
linetemp = lines[eachind+iatom].strip().split()
#print(linetemp)
linetemp = [eval(each) for each in linetemp]
for icolnm in range(Ncolnm):
normalmodes[ntriv+nmflag+icolnm][iatom] = linetemp[2+icolnm*3:5+icolnm*3]
nmflag += Ncolnm
return normalmodes, freqs, ele, linear
def nm2cart(nms,atom_mass):
Natoms = len(nms[0])
cart_nms = np.copy(nms)
for inm in range(len(nms)):
# reduced mass
reduced_amu = 0.0
for iatom in range(Natoms):
reduced_amu += np.sum(nms[inm][iatom]**2*atom_mass[iatom])
# cartesian normal modes (au)
for iatom in range(Natoms):
cart_nms[inm][iatom] /= np.sqrt(reduced_amu*1822.888515)
# Mass weighted normal modes:
w_nmode = np.copy(cart_nms)
for inm in range(len(nms)):
for iatom in range(Natoms):
w_nmode[inm][iatom] *= np.sqrt(atom_mass[iatom]*1822.888515)
return cart_nms, w_nmode
def herm2(xarg,x,nn,cost):
h = np.zeros(55)
h[0] = 1.0
h[1] = xarg + xarg
if (nn>=2):
for mm in range(1,nn+1):
ff = xarg*h[mm] - (mm-1)*h[mm-1]
h[mm+1] = ff+ff
ff = h[nn]
herm1 = ff*cost*np.exp(-x*x*0.5)
value = herm1*herm1
return value
# NewtonX-2.4-B06 initcond/sample.f90
def sample(amp,ivib,quant,shift,broad):
turnpo = np.sqrt(ivib*2.0+1.0)
xmin = -turnpo - quant
xdel = -xmin - xmin
# Initial conditions in the vibrational ground state
if ivib == 0:
while True:
x1 = np.random.random()
xx = xmin + xdel*x1
fsval = np.exp(-xx*xx)
x2 = np.random.random()
if (x2 < fsval):
break
xarg = xx / np.sqrt(broad) + shift
coor = xarg*amp
return coor
# Initial condition in vibrational excited states
nstep = 50*ivib
dlt = 1.0/nstep
fzb = -1000.0
k = 0
xarg = 0.0
while (xarg < -xmin):
k = k+1
xx = xmin+dlt*k
xarg = xx / np.sqrt(broad) + shift
fz = herm2(xarg,xx,ivib,1.0)
if (fz >= fzb): fzb = fz
fmax = fzb
cost = np.sqrt(1.0/fmax)
while True:
x1 = np.random.random()
xx = xmin + xdel*x1
xarg = xx / np.sqrt(broad) + shift
fsval = herm2(xarg,xx,ivib,cost)
x2 = np.random.random()
if (x2<= fsval): break
coor = xarg*amp
return coor
# NewtonX-2.4-B06 initcond/initqp.f90
# !! the len of all normal mode related array is numcoo
def rdmol(atom_mass,linear,freq):
eps = 1.0E-15
atmau = np.array(atom_mass)*1822.888515
Natoms = len(atom_mass)
numcoo = 3*Natoms - 6 + linear
anq = np.zeros(numcoo) # vibrational quantum numbers (all zero)
freq_au = np.array(freq) * 4.55633539E-6 # cm 2 au
# factor used for quantum distribution
amp = np.zeros(3*Natoms)
for i in range(numcoo):
amp[i] = 1.0 / np.sqrt(freq_au[i])
return anq, amp, freq_au, atmau
# NewtonX-2.4-B06 initcond/inqp.f90
def inqp(geomEq,cn,w,anq,amp,numcoo,atmau):
eps = 1.0E-15
cint = []
dcint = []
dum0 = 0.0
dum1 = 0.0
Natoms = len(atmau)
samp_points = np.zeros((2,numcoo))
q = np.array(geomEq) # size: (Natoms,3)
v = np.zeros((Natoms,3))
for i in range(numcoo):
wwai = amp[i]*w[i]
ivb = int(abs(anq[i])+0.5)
quant = 4.0
cint.append(sample(amp[i],ivb,quant,0.0,1.0))
dcint.append(sample(wwai,ivb,quant,0.0,1.0))
enmi = w[i]*(ivb+0.5) # Vibrational energies
epot = 0.5*(w[i]*cint[i])**2 # Potential energies
samp_points[0][i] = cint[i]/amp[i]
samp_points[1][i] = dcint[i]/wwai
dum0 += epot
dum1 += enmi
for iatom in range(Natoms):
for icoord in range(3):
fac = 1.0/np.sqrt(atmau[iatom])
for i in range(numcoo):
q[iatom][icoord] += cn[i][iatom][icoord]*cint[i]*fac
v[iatom][icoord] += cn[i][iatom][icoord]*dcint[i]*fac
return q,v # Unit: a.u.
def wignersample(npoints,molecule):
qlist = []
vlist = []
geomEq = np.array([each.xyz_coordinates for each in molecule.atoms])
mass = molecule.get_nuclear_masses()
mass = mass.reshape(1,len(mass))
# Calculate normal modes from Hessian matrix
#nm,freq,ele,linear_int = readGaussianNM(nmfile)
if not 'hessian' in molecule.__dict__.keys():
stopper.stopMLatom('Hessian matrix not found -- cannot do wigner sampling')
freq,nm,_,_ = thermo.vibrational_analysis(mass,molecule.hessian,mode_type='MDU')
linear = molecule.is_it_linear()
Natoms = len(molecule.atoms)
if linear:
ntriv = 5
linear_int = 1
else:
ntriv = 6
linear_int = 0
freq = freq[ntriv:]
for itriv in range(ntriv):
nm[itriv] = 1.0 / np.sqrt(3*Natoms)
numcoo = len(freq)
ntriv = len(nm)-numcoo
#print(len(nm),len(freq),len(ele),linear)
atom_mass = [each.nuclear_mass for each in molecule.atoms]
cart_nms, w_nmode = nm2cart(nm,atom_mass)
geomEq_au = np.array(geomEq) / constants.Bohr2Angstrom
anq, amp, freq_au, atmau = rdmol(atom_mass,linear_int,freq)
for ipoint in range(npoints):
q,v = inqp(geomEq_au,w_nmode[ntriv:],freq_au,anq,amp,numcoo,atmau)
qlist.append(q)
vlist.append(v)
# Transform from a.u. to Angstrom & fs
qlist = np.array(qlist) * constants.Bohr2Angstrom
vlist = np.array(vlist) * constants.Bohr2Angstrom / constants.au2fs
return qlist, vlist, linear_int
# generate random velocities without angular momentum for linear molecule
def generate_random_velocities_for_linear_molecule(molecule):
np.random.seed()
Natoms = len(molecule.atoms)
coord = molecule.xyz_coordinates
randnum = np.random.randn(Natoms)
avgnum = np.average(randnum)
randnum = randnum - avgnum
vec = coord[1] - coord[0]
rand_velocities = [vec*each for each in randnum]
return rand_velocities