#* **************************************************************************
#*
#* CDTK, Chemical Dynamics Toolkit
#* A modular system for chemical dynamics applications and more
#*
#* Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016
# H_1 ---- C_a ===== C_b ---- H_2
#*
#* Copyright (C) 2017, 2018, 2019
#* Ralph Welsch, DESY, <ralph.welsch@desy.de>
#*
#* Copyright (C) 2020, 2021, 2022, 2023
#* Ludger Inhester, DESY, ludger.inhester@cfel.de
#*
#* This file is part of CDTK.
#*
#* CDTK is free software: you can redistribute it and/or modify
#* it under the terms of the GNU General Public License as published by
#* the Free Software Foundation, either version 3 of the License, or
#* (at your option) any later version.
#*
#* This program is distributed in the hope that it will be useful,
#* but WITHOUT ANY WARRANTY; without even the implied warranty of
#* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#* GNU General Public License for more details.
#*
#* You should have received a copy of the GNU General Public License
#* along with this program. If not, see <http://www.gnu.org/licenses/>.
#*
#* **************************************************************************
import math as mt
import numpy as np
import CDTK.Tools.Conversion as cv
########################################################################################
#
# H_1 ---- C_a ===== C_b ---- H_2
#
# Description of the [C2H2]3+ system in terms of three singly charged fragments:
# CH+: H_1 ---- C_a
# C+: C_b
# H+: H_2
########################################################################################
# Parameters of the model
d = 3.7349*cv.ev2au # dissociation energy of CH+
a = 1.04488 # Curvature of Morse pot.
re = 1.136*cv.an2au # Equilibrium distance of Morse Pot.
q_a = 0.68 # partial charge of C_a
q_1 = 1.0 - q_a # partial charge of H_1
[docs]
def pot(X):
"""
X[0],X[1],X[2] C_a
X[3],X[4],X[5] C_b
X[6],X[7],X[8] H_1
X[9],X[10],X[11] H_2
"""
# interatomic distances
r1a = np.linalg.norm( X[6:9] - X[0:3] )
r1b = np.linalg.norm( X[6:9] - X[3:6] )
r12 = np.linalg.norm( X[6:9] - X[9:12] )
rab = np.linalg.norm( X[0:3] - X[3:6] )
ra2 = np.linalg.norm( X[0:3] - X[9:12] )
rb2 = np.linalg.norm( X[3:6] - X[9:12] )
v = 0.0
v = v + q_1/r1b + q_1/r12 + q_a/rab + q_a/ra2 + 1.0/rb2
v = v + _morse(r1a)
return v
[docs]
def grd(X):
"""
X[0],X[1],X[2] C_a
X[3],X[4],X[5] C_b
X[6],X[7],X[8] H_1
X[9],X[10],X[11] H_2
"""
# interatomic distances
r1a = np.linalg.norm( X[6:9] - X[0:3] )
r1b = np.linalg.norm( X[6:9] - X[3:6] )
r12 = np.linalg.norm( X[6:9] - X[9:12] )
rab = np.linalg.norm( X[0:3] - X[3:6] )
ra2 = np.linalg.norm( X[0:3] - X[9:12] )
rb2 = np.linalg.norm( X[3:6] - X[9:12] )
# interatomic vectors
R1a = X[6:9] - X[0:3]
R1b = X[6:9] - X[3:6]
R12 = X[6:9] - X[9:12]
Rab = X[0:3] - X[3:6]
Ra2 = X[0:3] - X[9:12]
Rb2 = X[3:6] - X[9:12]
G = np.zeros(12,float)
# Coulomb terms
# 1-b
f1b = q_1/r1b**3 * R1b
G[6:9] = G[6:9] - f1b
G[3:6] = G[3:6] + f1b
# 1-2
f12 = q_1/r12**3 * R12
G[6:9] = G[6:9] - f12
G[9:12] = G[9:12] + f12
# a-b
fab = q_a/rab**3 * Rab
G[0:3] = G[0:3] - fab
G[3:6] = G[3:6] + fab
# a-2
fa2 = q_a/ra2**3 * Ra2
G[0:3] = G[0:3] - fa2
G[9:12] = G[9:12] + fa2
# b-2
fb2 = 1.0/rb2**3 * Rb2
G[3:6] = G[3:6] - fb2
G[9:12] = G[9:12] + fb2
# Morse
e = mt.exp(-a*(r1a-re))
f1a = 2.0*d*a/r1a * (e-e**2) * R1a
G[6:9] = G[6:9] + f1a
G[0:3] = G[0:3] - f1a
return G
def _morse(r):
e = mt.exp(-a*(r-re))
return d*(1.0 - e)**2
def _numgrd(x):
# numerical gradient, for testing only
dx = 0.001 # step for numerical differentiation
grad = np.zeros(x.size,float)
for i in range(x.size):
x1 = x.copy()
x1[i] = x1[i] + 0.5*dx
f1 = pot(x1)
x0 = x.copy()
x0[i] = x0[i] - 0.5*dx
f0 = pot(x0)
grad[i] = (f1-f0)/dx
return grad
def _remap_coordinates(x0,v0=None):
"""
Swap hydrogens if necessary such that the closest to a C feels the Morse pot.
C_a and H_1 are those connected by a Morse potential.
x[0],x[1],x[2] C_a
x[3],x[4],x[5] C_b
x[6],x[7],x[8] H_1
x[9],x[10],x[11] H_2
"""
x0.shape = (4,3)
if v0 is not None:
v0.shape = (4,3)
isvelocity = True
else:
v0 = np.zeros((4,3),float)
isvelocity = False
r1a = np.linalg.norm( x0[2] - x0[0] )
r1b = np.linalg.norm( x0[2] - x0[1] )
ra2 = np.linalg.norm( x0[0] - x0[3] )
rb2 = np.linalg.norm( x0[1] - x0[3] )
x = x0.copy()
v = v0.copy()
rvec = np.array([r1a,r1b,ra2,rb2])
imin = rvec.argmin()
if imin == 1: # exchange carbons
x[0] = x0[1]
x[1] = x0[0]
v[0] = v0[1]
v[1] = v0[0]
if imin == 2: # exchange hydrogens
x[2] = x0[3]
x[3] = x0[2]
v[2] = v0[3]
v[3] = v0[2]
if imin == 3: # exchange carbons and hydrogens
x[0] = x0[1]
x[1] = x0[0]
x[2] = x0[3]
x[3] = x0[2]
v[0] = v0[1]
v[1] = v0[0]
v[2] = v0[3]
v[3] = v0[2]
x.shape = (12,)
v.shape = (12,)
if isvelocity:
return x,v
else:
return x