#* **************************************************************************
#*
#* CDTK, Chemical Dynamics Toolkit
#* A modular system for chemical dynamics applications and more
#*
#* Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016
#* Oriol Vendrell, DESY, <oriol.vendrell@desy.de>
#*
#* 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/>.
#*
#* **************************************************************************
from __future__ import unicode_literals
import sys
import os
import shutil
from optparse import OptionParser
import imp
import numpy as np
import CDTK.Dynamics.Trajectory_SH as trjsh
import CDTK.Tools.Inputfile as inf
import CDTK.Tools.Conversion as cv
import CDTK.Tools.Utils as uti
import CDTK.Models
[docs]
def terminate(m):
sys.stderr.write('\n It seems something went wrong in xmodel:\n')
sys.exit(m)
[docs]
def start():
"""
Wrapper for starting the program. Needed for good unit test coding.
"""
IDLEN = 8
TINY = 1.0e-8
PREFIX = 'run'
# --------------------------------------------------------------------------
# Parse command line options
# --------------------------------------------------------------------------
parser=OptionParser()
parser.add_option('-i','--input-file',
dest='input_file',
type='str',
default=None,
help='Path to input file')
parser.add_option('-d','--project-dir',
dest='project_dir',
type='str',
default=None,
help='Project name used to for output directory')
opts, args = parser.parse_args(sys.argv[1:])
# --------------------------------------------------------------------------
# Parse input file
# --------------------------------------------------------------------------
basename = opts.input_file.split('.')[0]
I = inf.Inputfile(opts.input_file).sections
# --------------------------------------------------------------------------
# Set projectname
# --------------------------------------------------------------------------
projdir = opts.project_dir
if not projdir:
_id = uti.getUniqueID(IDLEN)
_date = uti.getDateString()
projdir = PREFIX + '_' + basename + '_' + _date + '_' + _id
else:
projdir = PREFIX + '_' + projdir
# --------------------------------------------------------------------------
# Create name directory, copy input file and change to it
# --------------------------------------------------------------------------
if not os.path.exists(projdir):
os.mkdir(projdir)
shutil.copy(opts.input_file,projdir+'/input')
os.chdir(projdir)
open('stop','w').close() # empty "stop" file
# --------------------------------------------------------------------------
# Read positions and velocities
# --------------------------------------------------------------------------
X0 = I['inicond'].get('pos',None)
V0 = I['inicond'].get('vel',None)
if X0 is None:
terminate('initial positions not found')
else:
X0 = np.asarray(X0,float)
if V0 is None:
V0 = np.zeros(len(X0),float)
else:
V0 = np.asarray(V0,float)
ncoord = len(X0)
# --------------------------------------------------------------------------
# Load model module
# --------------------------------------------------------------------------
modelname = I['system']['model'][0]
f=imp.find_module(modelname,CDTK.Models.__path__)
model = imp.load_module(modelname,f[0],f[1],f[2])
# --------------------------------------------------------------------------
# Initialize trajectory object
# --------------------------------------------------------------------------
dt = cv.numval( I['trajectory']['dt'] , 'time' )
tf = cv.numval( I['trajectory']['tf'] , 'time' )
try:
rseed = int( I['system']['rseed'][0] )
except:
rseed = None
if I['quantum']['type'][0] == 'fssh':
nstates = int( I['quantum']['nstates'][0] )
istate = int( I['quantum']['istate'][0] )
trajectory = trjsh.Trajectory_SH(ncoord,nstates)
trajectory.rseed = rseed
trajectory.S = istate
trajectory.C[:] = 0
trajectory.C[istate]= 1
trajectory.DIR = os.getcwd()
trajectory.R = X0
trajectory.V = V0
trajectory.M = model.masses
trajectory.dt = dt
trajectory.f_EGrad = model.f_EGrad
trajectory.f_W = model.f_W
trajectory.f_D = model.f_D
rescaling = I['quantum'].get('rescaling',None)
if rescaling:
if rescaling[0] == 'nac':
trajectory.rescaling = 'nac'
trajectory.f_NAC = model.f_NAC
if rescaling[0] == 'grd':
trajectory.rescaling = 'grd'
else:
terminate('$quantum|type needs to be specified')
# --------------------------------------------------------------------------
# Run trajectory
# --------------------------------------------------------------------------
while trajectory.t < tf:
if os.path.isfile('stop'):
trajectory.integrate() # dt step
else:
break
os.chdir('..')
if __name__ == "__main__":
start()