#*********************************************************************************
#*
#* Copyright (C) 2020 Mohammad Shakiba, Brendan Smith, Alexey V. Akimov
#* This file is distributed 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.
#* See the file LICENSE in the root directory of this distribution
#* or <http://www.gnu.org/licenses/>.
#*
#*********************************************************************************/
"""
.. module:: step2_many_body
:platform: Unix, Windows
:synopsis: This module implements functions for computing the overlap and nonadiabatic couplings matrices.
.. moduleauthors::
Mohammad Shakiba, Brendan Smith, Alexey V. Akimov
"""
import numpy as np
import math
import sys
import os
import time
import copy
import multiprocessing as mp
# Concurrency for large systems
import concurrent.futures
from liblibra_core import *
from libra_py import data_conv
from libra_py import cube_file_methods
from libra_py import CP2K_methods
from libra_py import Gaussian_methods
from libra_py import DFTB_methods
from libra_py.workflows.nbra import mapping
from libra_py.workflows.nbra import step3
from libra_py import units
import util.libutil as comn
# This file is temp only. These functions will eventually be placed in Libra somewhere ...
[docs]def curr_and_final_step_job( istep, fstep, njobs, njob ):
"""
This function is used to determine the initial and final step of a job when distributing
the molecular dynamics trajectory over a number of jobs.
Args:
istep (int): The initial time step for molecular dynamics trajectory which
the user wants to start the calculations with. This paramter
is set in the input file for submitting the jobs.
fstep (int): The final time step for molecular dynamics trajectory which
the user wants to start the calculations with. This paramter
is set in the input file for submitting the jobs.
njobs (int): The number of jobs specified by user.
njob (int): The job number.
Returns:
job_init_step (int): The job initial time step.
job_final_step (int): The job final time step.
"""
total_steps = fstep - istep + 1
nsteps_per_job = int(total_steps/njobs)
if njob == 0:
job_init_step = istep
job_final_step = istep + nsteps_per_job - 1
elif njob > 0 and njob < ( njobs - 1 ):
job_init_step = istep + njob * nsteps_per_job - 1
job_final_step = istep + ( njob + 1 ) * nsteps_per_job - 1
elif njob == ( njobs - 1 ):
job_init_step = istep + njob * nsteps_per_job - 1
job_final_step = fstep
return job_init_step, job_final_step
[docs]def normalize_ci_coefficients(ci_coefficients_raw_unnorm):
"""
This funciton normalizes the list of configuration interaction (CI) coefficients.
Args:
ci_coefficients_raw_unnorm (list): The list containing the lists of unnormalized CI coefficients.
Returns:
ci_coefficients_raw_norm (list): The list containing the lists of normalized CI coefficients.
"""
# Number of states contributing in the excited state
nstates = len(ci_coefficients_raw_unnorm)
# Creating an empty list to store the normalized CI coefficients
ci_coefficients_raw_norm = []
for i in range(nstates):
#### ordinary way without using numpy
# Set up an initial parameter to compute the norm of the vector
norm = 0.0
# For each list of CI coefficients
ci_coefficients_raw_norm.append( [] )
for j in ci_coefficients_raw_unnorm[i]:
# sum of CI coefficients square
norm += j*j
# Compute the norm by taking the square root of the sum to compute the norm
norm = math.sqrt(norm)
# numpy way to compute the normalized CI coefficients
ci_coefficients_raw_unnorm[i] = np.array(ci_coefficients_raw_unnorm[i])
ci_coefficients_raw_norm[i] = ci_coefficients_raw_unnorm[i] / norm
ci_coefficients_raw_norm[i] = list(ci_coefficients_raw_norm[i])
return ci_coefficients_raw_norm
[docs]def get_excitation_analysis_output(params):
"""
This function reads the information of the excited states from the log file of the single point calculations.
Args:
params (dict):
logfile_directory (str): The log files directory.
es_software (str): The name of the software used to calculate the energy calculations.
curr_step (int): The current time step of the calculations.
isUKS (int): This parameter is set for spin restricted and unrestricted calculations. When it is
set to 1 it means that unrestricted calculations were set in the input file otherwise
it is restricted.
Returns:
excitation_energies (numpy array): The excitation energies of the curr_step.
ci_basis_raw (numpy array): The excited states which contains the occupied and virtual orbitals.
ci_coefficients_raw_unnorm (numpy array): The excited states CI coefficients.
spin_components (numpy array): Contains the excited states spin components ('alp' for alpha spin and 'bet' for beta spin)
"""
critical_params = [ "curr_step" ]
# Default parameters
default_params = { "isUKS": 0, "es_software": "cp2k", "logfile_directory": "logfiles" }
# Check input
comn.check_input(params, default_params, critical_params)
logfile_directory = params["logfile_directory"]
es_software = params["es_software"]
curr_step = params["curr_step"]
if es_software == "cp2k":
logfile_name = F"{logfile_directory}/step_{curr_step}.log"
params.update({"logfile_name":logfile_name})
# Extract the excitation energies and their configurations from the output file.
excitation_energies, ci_basis_raw, ci_coefficients_raw_unnorm, spin_components = \
CP2K_methods.read_cp2k_tddfpt_log_file( params )
elif es_software == "gaussian":
logfile_name = F"{logfile_directory}/step_{curr_step}.log"
params.update({"logfile_name":logfile_name})
excitation_energies, ci_basis_raw, ci_coefficients_raw_unnorm, spin_components = \
Gaussian_methods.read_gaussian_tddft_log_file( params )
elif es_software == "dftb+":
logfile_name = F"{logfile_directory}/TRA_{curr_step}.DAT"
params.update({"logfile_name":logfile_name})
excitation_energies, ci_basis_raw, ci_coefficients_raw_unnorm, spin_components = \
DFTB_methods.read_dftbplus_TRA_file( params )
return excitation_energies, ci_basis_raw, ci_coefficients_raw_unnorm, spin_components
[docs]def integrate_cube_set( cubefiles_set_1, cubefiles_set_2, dv ):
"""
This function comutes the overlap matrix between two set of cube files.
Args:
cubefiles_set_1 (list): The list of cube files of the curr_step.
cubefiles_set_2 (list): The list of cube files of the previous step.
dv (float): The integration element obtained from the cube files.
Returns:
overlap_matrix (numpy 2D array): The overlap between the two set of cube files.
"""
# Initialize the overlap matrix
# Note: The overlap matrix is a square matrix so the cubefiles_set_1 and cubefiles_set_2 must have the same length
overlap_matrix = np.zeros( ( len( cubefiles_set_2 ) , len( cubefiles_set_1 ) ) )
# The overlaps between the cube files of the first set and second set
for i in range( 0, len( cubefiles_set_2 ) ):
for j in range( 0, len( cubefiles_set_1 ) ):
# Use the integrate_cube to compute the integration between the wavefunctions.
overlap_matrix[i][j] = cube_file_methods.integrate_cube( cubefiles_set_2[i], cubefiles_set_1[j], dv )
return overlap_matrix
[docs]def compute_cube_ks_overlaps( cubefiles_prev, params):
"""
This function computes overlaps between cube files of two time steps. In order to not read the cube files
twice it returns the cube files of the current step, and gets the cube files of the previous step.
Args:
cubefiles_prev (list): The list containing th cube files of the previous step.
params (dict):
curr_step (int): The current time step.
isUKS (int): This parameter is set for spin restricted and unrestricted calculations. When it is
set to 1 it means that unrestricted calculations were set in the input file otherwise
it is restricted.
nprocs (int): The number of processors used to read the cube files and perform the integration.
Returns:
cubefiles_curr (list): The list of the read current step cube files.
S_ks_prev (2D numpy array): The overlap matrix of the wavefunctions for the previous time step.
S_ks_curr (2D numpy array): The overlap matrix of the wavefunctions for the current time step.
St_ks (2D numpy array): The overlap matrix between wavefunctions of the two time step.
"""
critical_params = [ "curr_step", "nprocs" ]
# Default parameters
default_params = { "isUKS": 0, "es_software": "cp2k"}
# Check input
comn.check_input(params, default_params, critical_params)
# Extract the variables
curr_step = int(params["curr_step"])
isUKS = int(params["isUKS"])
nprocs = int(params["nprocs"])
es_software = params["es_software"]
print('---------------------------------------------------------')
print('Starting the calculations by reading the cube files \n')
print(' Step time %i\n'%(curr_step))
print('-------------------------------------------------------\n')
print('Reading cube files with %i number of processors'%nprocs)
# First read all the .cube files (This is the most time consuming part)
# Creating a pool and assigning the number of processors for mp.Pool
pool = mp.Pool(processes=nprocs)
# for curr_step
print("Reading the cubes for step ", curr_step)
# generate the cube file names produced by CP2K
cubefile_names_curr = CP2K_methods.cube_file_names_cp2k( params )
#for cubefile in cubefile_names_curr:
# os.system( "/gpfs/scratch/brendan/cp2k/tools/cubecruncher/cubecruncher.x -center geo -i %s -o %s-1.cube " % ( cubefile, cubefile.replace( ".cube", "" ) ) )
# os.system( "rm %s" % cubefile)
# os.system( "mv %s-1.cube %s" % ( cubefile.replace(".cube",""), cubefile ) )
# Reading the cube files
cubefiles_curr = []
# Apply pool.map to the cube_file_methods to the set of variables of the cubefile_names_curr
cubefiles_curr = pool.map( cube_file_methods.read_cube, cubefile_names_curr )
# Close the pool
pool.close()
# Calculate the dv element for integration
dv = cube_file_methods.grid_volume( cubefile_names_curr[0] )
#### Compute the overlaps
# If spin unrestricted calculations were set to run
if isUKS == 1:
# The cube files for alpha and beta spin for the previous time step
# The alpha cubes are the even indices of the read cube files
alp_cubes_prev = cubefiles_prev[0::2]
# The beta cubes are the odd indices of the read cube files
bet_cubes_prev = cubefiles_prev[1::2]
# The cube files for alpha and beta spin for the current time step
# The alpha cubes are the even indices of the read cube files
alp_cubes_curr = cubefiles_curr[0::2]
# The beta cubes are the odd indices of the read cube files
bet_cubes_curr = cubefiles_curr[1::2]
# Initializing the overlap matrices for alpha and beta spin with zero matrices
zero_mat_alp = np.zeros( ( len( alp_cubes_prev ), len( alp_cubes_curr ) ) )
zero_mat_bet = np.zeros( ( len( bet_cubes_prev ), len( bet_cubes_curr ) ) )
### Using concurrency to compute the integration and overlap matrices
with concurrent.futures.ThreadPoolExecutor(max_workers=nprocs) as executor:
# <psi_alp(t-1)|psi_alp(t)>
int_1 = executor.submit(integrate_cube_set, alp_cubes_prev, alp_cubes_curr, dv )
# <psi_bet(t-1)|psi_bet(t)>
int_2 = executor.submit(integrate_cube_set, bet_cubes_prev, bet_cubes_curr, dv )
# <psi_alp(t-1)|psi_alp(t-1)>
int_3 = executor.submit(integrate_cube_set, alp_cubes_prev, alp_cubes_prev, dv )
# <psi_bet(t-1)|psi_bet(t-1)>
int_4 = executor.submit(integrate_cube_set, bet_cubes_prev, bet_cubes_prev, dv )
# <psi_alp(t)|psi_alp(t)>
int_5 = executor.submit(integrate_cube_set, alp_cubes_curr, alp_cubes_curr, dv )
# <psi_bet(t)|psi_bet(t)>
int_6 = executor.submit(integrate_cube_set, bet_cubes_curr, bet_cubes_curr, dv )
# Extracting the results for each of the submitted jobs
St_alp_alp = int_1.result()
St_bet_bet = int_2.result()
# The overlap between cube files at times t-1 and t
S_alp_alp_prev = int_3.result()
S_bet_bet_prev = int_4.result()
S_alp_alp_curr = int_5.result()
S_bet_bet_curr = int_6.result()
else:
### Using the concurrency
with concurrent.futures.ThreadPoolExecutor(max_workers=nprocs) as executor:
# <psi(t-1)|psi(t-1)>
int_1 = executor.submit(integrate_cube_set, cubefiles_prev, cubefiles_prev, dv )
# <psi(t)|psi(t)>
int_2 = executor.submit(integrate_cube_set, cubefiles_curr, cubefiles_curr, dv )
# <psi(t-1)|psi(t)>
int_3 = executor.submit(integrate_cube_set, cubefiles_prev, cubefiles_curr, dv )
# Extracting the results
S_prev = int_1.result()
S_curr = int_2.result()
St = int_3.result()
# These are used to form the block matrices for two-spinor format
zero_mat_alp = np.zeros( ( len( S_prev ), len( S_curr ) ) )
zero_mat_bet = np.zeros( ( len( S_prev ), len( S_curr ) ) )
# The alpha and beta spin have the same overlap matrix in spin restricted case
S_alp_alp_prev, S_bet_bet_prev = S_prev, S_prev
S_alp_alp_curr, S_bet_bet_curr = S_curr, S_curr
St_alp_alp, St_bet_bet = St, St
# Storing the overlap matices in two-spinor format by forming the block matrices
S_ks_prev = data_conv.form_block_matrix( S_alp_alp_prev, zero_mat_alp, zero_mat_bet, S_bet_bet_prev )
S_ks_curr = data_conv.form_block_matrix( S_alp_alp_curr, zero_mat_alp, zero_mat_bet, S_bet_bet_curr )
St_ks = data_conv.form_block_matrix( St_alp_alp, zero_mat_alp, zero_mat_bet, St_bet_bet )
return cubefiles_curr, S_ks_prev, S_ks_curr, St_ks
[docs]def reindex_cp2k_sd_states( ks_orbital_homo_index, ks_orbital_indicies, sd_basis_states, sd_format=2 ):
"""
ks_orbital_homo_index: Index of the homo ks orbital, from 1
ks_orbital_indicies: Range of the considered ks orbtials. Ex) [8,9,10,11], where 9 is homo orbtial index (from 1)
sd_basis_states( list of lists of lists ): A list of Slater determinants, where each slater determinant is a excitation in the Kohn-Sham
basis. This function assumes that all Kohn-Sham excitations are for alpha electrons. To
differentiate between alpha and beta excitations, elements of sd_basis_states contain spin
information.
Ex) sd_basis_states[0] = [ [9,10], "alp" ]
sd_basis_states[1] = [ [9,10], "bet" ]
Returns:
Reindexed Slater determinant basis states in terms of the ks_orbital_homo_index. The SD states returned can be either in one of two ways
1. All alpha orbitals before beta orbtials. This format is still a fixed slot format. It is used when used both alpha and beta spin channels
Ex.) ks_orbital_homo_index = 68, ks_orbital_indicies = [67, 68, 69, 70, 71, 72]
The ground state will be: [1, 2, -7, -8]
[ [68,69], "alp" ] --> [1, 3, -7, -8]
[ [68,69], "bet" ] --> [1, 2, -7, -9]
2. Fixed slot format. this format used only the alpha spin channel
Ex.) ks_orbital_homo_index = 68, ks_orbital_indicies = [67, 68, 69, 70, 71, 72]
The ground state will be: [1, -1, 2, -2]
[ [68,69], "alp" ] --> [1, -1, 3, -2]
[ [68,69], "bet" ] --> [1, -1, 2, -3]
"""
print("\nEntered reindex_cp2k_sd_states function")
print("ks_orbital_homo_index", ks_orbital_homo_index)
print("ks_orbital_indicies", ks_orbital_indicies)
print("sd_basis_states", sd_basis_states)
# We need to update the indexing of the sd_basis - in terms of the rows and cols of St_KS
# reindex ks orbs according to the matrix size
n_alp_ks_orbs = len(ks_orbital_indicies)
alp_homo_matrix_index = 0
for i in range( n_alp_ks_orbs ):
if ks_orbital_indicies[i] == ks_orbital_homo_index:
alp_homo_matrix_index = i
print("alp_homo_matrix_index",alp_homo_matrix_index)
ks_orbs_new_index = []
for i in range( n_alp_ks_orbs ):
ks_orbs_new_index.append( i+1 )
# Form excited state SDs
excitations = []
# For each Slater determinant basis state, which could have spin_component "alp" or "bet"
for j in range( len( sd_basis_states ) ):
print(int( sd_basis_states[j][0][0] ))
print(int( sd_basis_states[j][0][1] ))
if sd_format == 1:
if sd_basis_states[j][1] == "alp":
initial_ks_orb = int( sd_basis_states[j][0][0] ) - ks_orbital_homo_index + alp_homo_matrix_index + 1
final_ks_orb = int( sd_basis_states[j][0][1] ) - ks_orbital_homo_index + alp_homo_matrix_index + 1
elif sd_basis_states[j][1] == "bet":
initial_ks_orb = int( sd_basis_states[j][0][0] ) - ks_orbital_homo_index + alp_homo_matrix_index + n_alp_ks_orbs + 1
final_ks_orb = int( sd_basis_states[j][0][1] ) - ks_orbital_homo_index + alp_homo_matrix_index + n_alp_ks_orbs + 1
elif sd_format == 2:
initial_ks_orb = int( sd_basis_states[j][0][0] ) - ks_orbital_homo_index + alp_homo_matrix_index + 1
final_ks_orb = int( sd_basis_states[j][0][1] ) - ks_orbital_homo_index + alp_homo_matrix_index + 1
if sd_basis_states[j][1] == "alp":
excitations.append( [initial_ks_orb, final_ks_orb] )
elif sd_basis_states[j][1] == "bet":
excitations.append( [-initial_ks_orb, -final_ks_orb] )
print( "excitations = ", excitations )
#sys.exit(0)
# Form ground-state SD first
sd_basis = [ [] ]
if sd_format == 1:
for i in range( 1, 2*len( ks_orbital_indicies ) ):
if i < alp_homo_matrix_index + 2:
sd_basis[0].append( i )
if i > n_alp_ks_orbs and i < n_alp_ks_orbs + alp_homo_matrix_index + 2:
sd_basis[0].append( -i )
elif sd_format == 2:
for i in range( 1, len( ks_orbital_indicies ) ):
if i < alp_homo_matrix_index + 2:
sd_basis[0].append( i )
sd_basis[0].append( -i )
print ( "ground state = ", sd_basis )
# Now that we have done the ground state slater
for j in range( len( excitations ) ):
sd_excitation = []
for sd_state in sd_basis[0]:
if sd_state==excitations[j][0]:
sd_excitation.append(excitations[j][1])
else:
sd_excitation.append(sd_state)
sd_basis.append( sd_excitation )
print ( sd_basis )
return sd_basis
[docs]def run_step2_many_body( params ):
"""
This function is the main function which runs the following calculations:
* Executes calls to perform CP2K, DFTB+, or Gaussian electronic structure calculations
* Computing the overlap, time-overlap, and energy eigenvalue matrices for the KS basis
* Molecular orbital visualization using VMD
Args:
params (dictionary):
min_band (integer): The minimum state number.
max_band (integer): The maximum state number.
ks_orbital_homo_index (integer): The index of the Kohn-Sham HOMO orbital (State numbers start from 1).
nsteps_this_job (integer): The number of steps in this job.
es_software (string): The electronic structure calculation software.
* Note: The current code works with CP2K but it can be generalized to other electronic structure
calculation software packages which are capable of producing cube files and performing
TD-DFT calculations such as Gaussian. So this parameter is for future implementations.
es_software_exe (string): The CP2K executable file or path.
es_software_input_template (str): The full path to CP2K input template.
project_name (string): The project name for CP2K input.
isUKS (integer): This is the flag for unrestricted spin calculations. When it is set to
1 it will perform the nrestricted calculations otherwise it will consider the restricted case.
trajectory_xyz_filename (string): The full path to trajectory xyz file.
njob (integer): The current job number.
nprocs (integer): The number of processors to be used.
res_dir (string): The full path to res directory where the output data are stored.
logfile_directory (string): The path to where the log files are stored.
Returns:
None
"""
# Current working directory
pwd = os.getcwd()
# Critical variables
critical_params = []
# Default parameters
default_params = { "min_band":1, "max_band":1, "ks_orbital_homo_index":0, "nsteps_this_job":1, 'trajectory_xyz_filename':"md.xyz", "isUKS": 0, "es_software": "cp2k", "es_software_exe": "cp2k.popt", "es_software_input_template": "cp2k_input_template.inp", "project_name": "Libra_CP2K", "njob": 1, "nprocs": 2, "logfile_directory": "logfiles", "istep": 0, "do_cube_visualization": 0, "states_to_be_plotted": [], "waveplot_exe":"/util/academic/dftbplus/20.2.1-arpack/bin/waveplot" }
# Check input
comn.check_input(params, default_params, critical_params)
# Extracting the min_band and max_band, we have to use int(...) since
# the numbers which are read from the bash are strings.
min_band = int( params["min_band"] )
max_band = int( params["max_band"] )
# Update params with min_band and max_band
params.update({"min_band":min_band})
params.update({"max_band":max_band})
# Define ks_orbital_indicies based on the min_band and max_band
ks_orbital_indicies = list( range( min_band, max_band+1 ) )
# Update params with ks_orbital_indicies
params.update({"ks_orbital_indicies":ks_orbital_indicies})
# The HOMO orbital index
ks_orbital_homo_index = int(params["ks_orbital_homo_index"])
# Number of steps for this job
nsteps_this_job = int(params["nsteps_this_job"])
es_software = params["es_software"].lower()
params.update( {"es_software": es_software} )
# The electronic structure software, for example cp2k
if es_software == "cp2k":
# The path to executable CP2K
cp2k_exe = params["es_software_exe"]
cp2k_input_template = params["es_software_input_template"]
elif es_software == "gaussian":
gaussian_exe = params["es_software_exe"]
gaussian_input_template = params["es_software_input_template"]
elif es_software == "dftb+":
dftbp_exe = params["es_software_exe"]
dftb_input_template = params["es_software_input_template"]
waveplot_exe = params["waveplot_exe"]
# The project name, this is important since the cube files are produced based on the project_name
project_name = params["project_name"]
# The path to the pre-computed trajectory
trajectory_xyz_filename = params["trajectory_xyz_filename"]
# The job number
njob = int(params["njob"])
# The number of processors
nprocs = int(params["nprocs"])
# Update params with nprocs
params.update({"nprocs":nprocs})
# The path to res directory where the overlap matrices, energies, and NACs are stored
res_dir = params["res_dir"]
# The path to log files produced by CP2K
logfile_directory = params["logfile_directory"]
# The current time step which here is set to istep
curr_step = int(params["istep"])
# Update the curr_step in params
params.update({ "curr_step": curr_step })
# Cube file visualization flag
do_cube_visualization = int(params["do_cube_visualization"])
# Make a directory for this job folder for storing the logfiles, cubefiles, and pdosfiles
os.mkdir("logfiles")
os.mkdir("cubefiles")
os.mkdir("pdosfiles")
if not os.path.exists("../../all_logfiles"):
os.mkdir("../../all_logfiles")
if not os.path.exists("../../all_pdosfiles"):
os.mkdir("../../all_pdosfiles")
#####################################################################################################
######################## Initializing lists for storing the data for each job #######################
#####################################################################################################
# For empty lists for the KS, SD, and CI bases. These are to be used for collecting data at each step
# Overlap of the wavefunctions at the same time step in this job
S_ks_job = []
# Overlap of the wavefunctions for two consecutive time steps in this job
St_ks_job = []
# Energy levels of the molecular orbitals for each time step in this job
E_ks_job = []
# The total energies for each time step in this job
total_energies_job = []
# The overlap matrices for Slater determinant basis for each time step in this job
S_sd_job = []
St_sd_job = []
####### Setting up the calculations for the initial step
timer1 = time.time()
# Set up an initial parameter to compute the norm of the vector
if es_software == "cp2k":
# Set up the cp2k input template with the atomic positions for the first timestep of this job batch
CP2K_methods.CP2K_input_static( cp2k_input_template, project_name, trajectory_xyz_filename, curr_step )
# Running the CP2K calculations
os.system("mpirun -np %d %s -i %s-%i.inp -o logfiles/step_%d.log "%( nprocs, cp2k_exe, project_name, curr_step, curr_step ) )
# After finishing the CP2K calculations move all the pdos and cube files to specified folders
os.system("mv *.pdos pdosfiles")
elif es_software == "gaussian":
Gaussian_methods.gaussian_input( project_name, curr_step, gaussian_input_template, trajectory_xyz_filename )
os.system("%s < %s-%i.gjf > logfiles/step_%d.log "%( gaussian_exe, project_name, curr_step, curr_step ) )
Gaussian_methods.cube_generator_gaussian(project_name,curr_step,ks_orbital_indicies[0],ks_orbital_indicies[-1],nprocs,'../../sample_cube_file.cube',int(params["isUKS"]))
elif es_software == "dftb+":
print( 'params["curr_step"]', params["curr_step"])
# Set up the dftb+ input template.
DFTB_methods.make_dftb_input( dftb_input_template, 0, params["curr_step"] )
# Convert the xyz file for currstep to the .gen file format
os.system("xyz2gen coord-"+str(params["curr_step"])+".xyz")
# Run the dftb+ program
os.system("srun %s"%( dftbp_exe ) )
os.system("mv band.out step_"+str(curr_step)+"_ks.log")
os.system("mv step_"+str(curr_step)+"_ks.log logfiles/.")
os.system("mv EXC.DAT EXC_"+str(curr_step)+".DAT")
os.system("mv EXC_"+str(curr_step)+".DAT logfiles/.")
os.system("mv TRA.DAT TRA_"+str(curr_step)+".DAT")
os.system("mv TRA_"+str(curr_step)+".DAT logfiles/.")
DFTB_methods.cube_generator_dftbplus( project_name, curr_step, ks_orbital_indicies[0], ks_orbital_indicies[-1], waveplot_exe, int(params["isUKS"]) )
os.system("mv *.cube cubefiles")
# Print the elapsed time for CP2K calculations for this step.
print(params["es_software"]," calculation time for step ", params["curr_step"]," was ", time.time() - timer1)
# We have compted the first SCF calculation for this job, now to read the output data and cubes
print("Reading the initial step using pool")
print("nprocs", nprocs)
# Creating a pool with nprocs
pool = mp.Pool( processes = nprocs )
# The cube file names produced by CP2K, Here we set it as prev since we
# don't want to read the cubes twice.
cubefile_names_prev = CP2K_methods.cube_file_names_cp2k( params )
print("Initial step for this job, the cube file names are:\n", cubefile_names_prev)
# We may have to "crunch" the cubes - this may be especially needed for periodic systems.
#for cubefile in cubefile_names_prev:
# os.system( "/gpfs/scratch/brendan/cp2k/tools/cubecruncher/cubecruncher.x -center geo -i %s -o %s-1.cube " % ( cubefile, cubefile.replace( ".cube", "" ) ) )
# os.system( "rm %s" % cubefile)
# os.system( "mv %s-1.cube %s" % ( cubefile.replace(".cube",""), cubefile ) )
cubefiles_prev = pool.map( cube_file_methods.read_cube, cubefile_names_prev )
# Call the pool to read the cube files
# Close the pool
pool.close()
# If we are to visualize the cubefiles, then we should compute the grid-mesh before hand to save time later
if do_cube_visualization == 1:
# This phase_factor_visual contains the phase factor correction for each state
# and is used to plot the cubes phase corrected
phase_factor_visual = np.ones( ( max_band - min_band + 1 ) * 2 )
# Update the params with phase_factor_visual
params.update({"phase_factor_visual":phase_factor_visual})
# The states which were set to be plotted
states_to_be_plotted = []
for state in list( params["states_to_be_plotted"].split(",") ):
states_to_be_plotted.append( int( state ) )
# Update the params for states_to_be_plotted
params.update({"states_to_be_plotted":states_to_be_plotted})
# Plot the cubes using VMD
cube_file_methods.plot_cubes( params )
# After reading the cubefiles for curr_step, we should delete them to be memory efficient
print("Removing cubefiles")
os.system("rm cubefiles/*")
# Froming the hvib_ks_re with energies for the first step in the job.
# This will extract the Kohn-Sham energies and total energy from the CP2K log files and forms the Hvib_real
hvib_ks_re, total_energy = form_Hvib_real( params )
# Appending the Kohn-Sham and total energies of this step in the E_ks_job and total_energies_job
E_ks_job.append( hvib_ks_re )
total_energies_job.append( total_energy )
curr_step += 1
#########################################
# All other steps after initial step for this job
for step in range( nsteps_this_job-1 ):
# Update the params with curr_step
params.update({ "curr_step":curr_step })
# Setting up the timer
timer1 = time.time()
if es_software == "cp2k":
# Creating the CP2K input
CP2K_methods.CP2K_input_static( cp2k_input_template, project_name, trajectory_xyz_filename, curr_step )
# Run CP2K
os.system("mpirun -np %d %s -i %s-%i.inp -o logfiles/step_%d.log "%( nprocs, cp2k_exe, project_name, curr_step, curr_step ) )
# Move all the pdos and cube files to one folder
os.system("mv *.pdos pdosfiles")
elif es_software == "gaussian":
Gaussian_methods.gaussian_input( project_name, curr_step, gaussian_input_template, trajectory_xyz_filename )
os.system("%s < %s-%i.gjf > logfiles/step_%d.log "%( gaussian_exe, project_name, curr_step, curr_step ) )
Gaussian_methods.cube_generator_gaussian(project_name,curr_step,ks_orbital_indicies[0],ks_orbital_indicies[-1],nprocs,'../../sample_cube_file.cube',int(params["isUKS"]))
elif es_software == "dftb+":
print( 'params["curr_step"]', params["curr_step"])
# Set up the dftb+ input template.
DFTB_methods.make_dftb_input( dftb_input_template, 0, params["curr_step"] )
# Convert the xyz file for currstep to the .gen file format
os.system("xyz2gen coord-"+str(params["curr_step"])+".xyz")
# Run the dftb+ program
os.system("srun %s"%( dftbp_exe ) )
os.system("mv band.out step_"+str(curr_step)+"_ks.log")
os.system("mv step_"+str(curr_step)+"_ks.log logfiles/.")
os.system("mv EXC.DAT EXC_"+str(curr_step)+".DAT")
os.system("mv EXC_"+str(curr_step)+".DAT logfiles/.")
os.system("mv TRA.DAT TRA_"+str(curr_step)+".DAT")
os.system("mv TRA_"+str(curr_step)+".DAT logfiles/.")
DFTB_methods.cube_generator_dftbplus( project_name, curr_step, ks_orbital_indicies[0], ks_orbital_indicies[-1], waveplot_exe, int(params["isUKS"]) )
os.system("mv *.cube cubefiles")
# Print the Timing
print("Elapsed time for step ", params["curr_step"]," was ", time.time() - timer1)
# Forming the hvib_ks_re the same as above
hvib_ks_re, total_energy = form_Hvib_real( params )
E_ks_job.append( hvib_ks_re )
total_energies_job.append( total_energy )
#============================== Computation of the overlap matrices=============================
# Now, read in the cube files for steps step-1 and step
# and form the S_ks and St_ks overlap matricies
cubefiles_curr, S_ks_prev, S_ks_curr, St_ks = compute_cube_ks_overlaps( cubefiles_prev, params )
# replace the cubefiles_prev with the cubefiles_curr
# with this we do not need to read the cubes twice
cubefiles_prev = cubefiles_curr
#===============================================================================================
# Plot the cube files using VMD
if do_cube_visualization == 1:
for row_index in range(len(St_ks)):
# Multiplying the phase factor for each states which needs to be plotted
if St_ks[row_index][row_index]<0:
phase_factor_visual[row_index] = phase_factor_visual[row_index] * (-1)
# Update the params with phase_factor_visual
params.update({"phase_factor_visual":phase_factor_visual})
# Plot the cubes
cube_file_methods.plot_cubes( params )
# Now remove all the cube files to be memory efficient
os.system("rm cubefiles/*")
# Appending the normalized overlap matrices in the _job variables
S_ks_job.append(S_ks_prev)
St_ks_job.append(St_ks)
if step == nsteps_this_job-2:
S_ks_job.append(S_ks_curr)
#np.savetxt(filename, np.array, fmt='%.16e', delimiter=" ")
np.savetxt("%s/S_ks_%d_re" % (res_dir, int(params["istep"])+step), S_ks_job[step], fmt='%.16e', delimiter=" ")
np.savetxt("%s/E_ks_%d_re" % (res_dir, int(params["istep"])+step), E_ks_job[step], fmt='%.16e', delimiter=" ")
np.savetxt("%s/St_ks_%d_re" % (res_dir, int(params["istep"])+step), St_ks_job[step], fmt='%.16e', delimiter=" ")
curr_step += 1
# Print out the KS Overlap and Energy matricies for the last step in this job batch
np.savetxt("%s/S_ks_%d_re" % (res_dir, int(params["istep"])+step+1), S_ks_job[step+1], fmt='%.16e', delimiter=" ")
np.savetxt("%s/E_ks_%d_re" % (res_dir, int(params["istep"])+step+1), E_ks_job[step+1], fmt='%.16e', delimiter=" ")
print("All steps were done successfully for this job!")
os.system("mv logfiles/* ../../all_logfiles/.")
os.system("mv pdosfiles/* ../../all_pdosfiles/.")
os.system("rm "+trajectory_xyz_filename)
os.system("rm *wfn")