nao.examples.qmd_C60 — NAO Examples: QMD C60#
A demonstration of the Nao module used to calculate the polarizability of C60 at rooms temperature using Quantum Molecular Dynamic (QMD). The relaxation and the ground state calculations are done using the Siesta program. This example will calculate the geometries of 5000 steps of C60 using Siesta and the polarizability will be obtained by calculating 200 steps. This example MUST be done on a server because it will perform 200 calculations. This example has been done for the Torque resources manager and will need to be adapted to your case before to run.
To start the calculations, we first need to generate the 5000 geometries using Molecular Dynamic (MD) with temperature controlled by means of a Nosé thermostat. For this we use the following siesta input file,
SystemLabel siesta
PAO.EnergyShift 50 meV
PAO.BasisSize DZP
PAO.SoftDefault .true.
PAO.BasisType split
XC.functional GGA
XC.authors PBE
WriteIonPlotFiles .false.
WriteCoorXmol .true.
WriteMDXmol .true.
%include "C60_org.fdf"
AtomCoorFormatOut Ang
DM.UseSaveDM .true.
DM.Tolerance 0.0001
DM.MixingWeight 0.25
MaxSCFIterations 400
DM.NumberPulay 4
MD.TypeOfRun Nose
MD.InitialTemperature 300 K
MD.TargetTemperature 300 K
MD.NumCGsteps 5000
MD.FinalTimeStep 5000
UseSaveData .true.
MeshCutOff 250 Ry
MD.UseSaveXV .true.
SaveRho .false.
COOP.Write .false.
WriteDenchar .false.
The file containing the position of the atoms can be downloaded with the following link.
- You run this calculation with the normal siesta command (don’t forget the C.psf file)::
siesta < siesta_C60_thermal.fdf > siesta.out
or with a bash script if you run on a server. Once, this first calculation finish you will be left with a siesta.ANI file containing the geometry at each time step. We will first generate the xyz files corresponding to the geometries for comomdity,
#!/bin/bash
fname=${1:-siesta.ANI}
suffi=${2:-0}
echo ".ANI file with MM geometries : " $fname
echo "Suffix to choose geometries for QM part : " $suffi
mkdir -p xyz
echo "cp $fname xyz/"
cp $fname xyz/
echo "cd xyz"
cd xyz
LIT=`head -n1 $fname`
NLINES_PER_FILE=$(expr $LIT + 2)
echo "rm x*"
rm x*
echo "split -l $NLINES_PER_FILE -d -a 6 $fname"
split -l $NLINES_PER_FILE -d -a 6 $fname
echo "xyz2fdf.py x0*$suffi > xyz2fdf.py.out"
xyz2fdf.py x0*$suffi &> xyz2fdf.py.out
cd ..
ls xyz/*.fdf
This will create a folder xyz containing the 5000 generated geometries. We need then to generate the fdf geometries file for siesta:
cd cyz
python xyz2fdf.py
xyz2fdf.py is a small python script that will generate the geometries for 200 calculations from the 5000 previous geometries because 5000 is a bit too much. The script look like,
from __future__ import division
import numpy as np
import ase.io as io
import subprocess
def write_geofdf(atoms, fname="geo.fdf"):
from ase.data import chemical_symbols
f=open(fname, "w")
f.write("NumberOfAtoms {0}\n".format(len(atoms)))
f.write("NumberOfSpecies {0}\n".format(len(set(atoms.numbers))))
f.write("\n")
f.write('%block ChemicalSpeciesLabel\n')
species_label = {}
for i, nb in enumerate(set(atoms.numbers)):
species_label[chemical_symbols[nb]] = i+1
f.write(' {0} {1} '.format(i+1, nb)+ chemical_symbols[nb]+'\n')
f.write("%endblock ChemicalSpeciesLabel\n")
f.write("\n")
f.write("AtomicCoordinatesFormat Ang\n")
f.write("%block AtomicCoordinatesAndAtomicSpecies\n")
for ia, atm in enumerate(atoms):
pos = atm.position
f.write("{0:.6f} {1:.6f} {2:.6f} {3} {4} ".format(pos[0], pos[1], pos[2],
species_label[atm.symbol], ia+1) + atm.symbol + "\n")
f.write("%endblock AtomicCoordinatesAndAtomicSpecies")
f.close()
xyz_range = np.arange(0, 5000, 25)
for i, xyz in enumerate(xyz_range):
if xyz < 10:
num = "00000{0}".format(xyz)
elif xyz < 100:
num = "0000{0}".format(xyz)
elif xyz < 1000:
num = "000{0}".format(xyz)
elif xyz < 10000:
num = "00{0}".format(xyz)
else:
raise ValueError("xyz too large?? {0}".format(xyz))
path = "calc_" + num
atoms = io.read("x"+num, format="xyz")
subprocess.call("mkdir " + path, shell=True)
write_geofdf(atoms, fname=path+"/geo.fdf")
subprocess.call("cp siesta_C60.fdf " + path, shell=True)
subprocess.call("cp C.psf " + path, shell=True)
The script will generate the siesta input files in different trajectory, don’t forget to add in the xyz directory the siesta input
SystemLabel siesta
PAO.EnergyShift 50 meV
PAO.BasisSize DZP
PAO.SoftDefault .true.
PAO.BasisType split
XC.functional GGA
XC.authors PBE
%include "geo.fdf"
AtomCoorFormatOut Ang
DM.MixingWeight 0.01
DM.Tolerance 0.0001
DM.MixingWeight 0.25
DM.NumberPulay 4
MaxSCFIterations 400
MD.TypeOfRun CG
MD.NumCGsteps 0
MD.MaxForceTol 0.02 eV/Ang
SolutionMethod Diagon
MeshCutOff 250 Ry
WriteCoorXmol .True.
WriteDenchar .True.
COOP.Write .True.
XML.Write .True.
- After running the xyz2fdf.py script you should have 200 folders containing the 3 input files,
siesta_C60.fdf
geo.fdf
C.psf
I advise you to create a new folder, and to copy all the calculations input into this new folder:
cd ..
mkdir calculations
cp -r xyz/calc_* calculations
cd calculations
Now we will need script to perform the siesta and nao calculations. Here, I will perform in total 50 jobs, each running 4 calculations. Each jobs will be using 6 cpus. I’m assuming that you are using the Torque managing system. We will use two python script to perform the calculations,
The first one is the script that you call in order to submit the 200 jobs and it call the second script calc_polarizability.py. This second script will perform the siesta and nao calculations for each configurations.
submit_jobs.py:
from __future__ import division
import numpy as np
import subprocess
# The bash script for Torque
script = """#!/bin/bash
#PBS -q parallel
#PBS -l nodes=1:ppn=6
#PBS -l mem=5gb
#PBS -l cput=400:00:00
#PBS -N calc_C60
ulimit -s unlimited
export NPROCS=`wc -l < $PBS_NODEFILE`
export OMP_NUM_THREADS=${NPROCS}
export PATH="/Path/To/Python/binary:$PATH"
export LD_LIBRARY_PATH="/Path/To/Needed/Library:$LD_LIBRARY_PATH" # libxc, libcint, libxcfun
export PYTHONPATH="/Path/To/Pyscf:$PYTHONPATH"
# ASE necessay for using ase.units
#ASE
export ASE_HOME=/Path/To/ASE
export PYTHONPATH="${ASE_HOME}:$PYTHONPATH"
export PATH="${ASE_HOME}/tools:$PATH"
# you may need to change the job directory
export LSCRATCH_DIR=/scratch/$USER/jobs/$PBS_JOBID
mkdir -p $LSCRATCH_DIR
cd $PBS_O_WORKDIR
# load the right module depending your pyscfi/siesta compilation
# you better to use the same compiler for both program to avoid problems
ml purge
ml load intel/2015b FFTW/3.3.4-intel-2015b
"""
# the range of our 200 calculations
xyz_range = np.arange(0, 5000, 25)
start = 0
step = 4
end = 0
while end < xyz_range.size:
end = start + step
if end > xyz_range.size:
end = xyz_range.size + 1
calcs = xyz_range[start:end]
include = ["calc_polarizability.py"]
for i, xyz in enumerate(calcs[0:calcs.shape[0]]):
if xyz < 10:
num = "00000{0}".format(xyz)
elif xyz < 100:
num = "0000{0}".format(xyz)
elif xyz < 1000:
num = "000{0}".format(xyz)
elif xyz < 10000:
num = "00{0}".format(xyz)
else:
raise ValueError("xyz too large?? {0}".format(xyz))
include.append("calc_"+num)
fcalc = calcs[0]
ecalc = calcs[calcs.shape[0]-1]+25
lines = script
for files in include:
lines += "cp -r " + files + " $LSCRATCH_DIR\n"
lines += "cd $LSCRATCH_DIR\n"
lines += "./calc_polarizability.py --np ${NPROCS} " + "--start {0} --end {1} >& calc_{0}to{1}.out\n".format(fcalc, ecalc)
lines += "export RESULTS_DIR=$PBS_O_WORKDIR\n"
lines += "mkdir -p $RESULTS_DIR\n"
lines += "cp -r * $RESULTS_DIR\n"
fname = "run.calc_C60_{0}to_{1}.sh".format(fcalc, ecalc)
f = open(fname, "w")
f.write(lines) # write bash script
f.close()
start = end
# submit job to Torque
subprocess.call("qsub " + fname, shell=True)
calc_polarizability.py
#!/Path/To/Python/binary
from __future__ import print_function, division
import os,numpy as np
from pyscf.nao import tddft_iter
import sys
import argparse
from ase.units import Ry, eV, Ha
import subprocess
from timeit import default_timer as timer
def run_tddft_iter():
t1 = timer()
td = tddft_iter(label="siesta", iter_broadening=0.0035/Ha, xc_code='LDA,PZ', level=0, tddft_iter_tol=1e-3, tol_loc=1e-4, tol_biloc=1e-6)
t2 = timer()
print("time tddft_iter = ", t2-t1)
omegas = np.arange(2.0, 8.0, 2.5E-3)/Ha + 1j*td.eps
t3 = timer()
pxx = -td.comp_polariz_nonin_ave(omegas)
t4 = timer()
print("time chi0 = ", t4-t3)
data = np.zeros((omegas.shape[0], 3))
data[:, 0] = omegas.real*Ha
data[:, 1] = pxx.real
data[:, 2] = pxx.imag
np.savetxt('polarizability_nonin_siesta.avg.txt', data)
t5 = timer()
pxx = -td.comp_polariz_inter_ave(omegas)
t6 = timer()
print("time chi = ", t6-t5)
data = np.zeros((omegas.shape[0], 3))
data[:, 0] = omegas.real*Ha
data[:, 1] = pxx.real
data[:, 2] = pxx.imag
np.savetxt('polarizability_inter_siesta.avg.txt', data)
print("nb call:")
print("rf0_ncalls = {0}, matvec_ncalls ={1}".format(td.rf0_ncalls, td.matvec_ncalls))
t7 = timer()
print("total time = ", t7-t1)
parser = argparse.ArgumentParser()
parser.add_argument('--np', type=int, default=1, help="number of processor to use")
parser.add_argument('--start', type=int, default=0, help="starting calc")
parser.add_argument('--end', type=int, default=25, help="end calculation")
args_par = parser.parse_args()
xyz_range = np.arange(args_par.start, args_par.end, 25)
mpi_exe="/scratch/mbarbry/intel/intelpython2/bin/mpirun"
siesta_path="/scratch/software/SIESTA/4.0b-485-intel-2015b/siesta"
siesta_exe = mpi_exe + " -np {0} ".format(args_par.np) + siesta_path + " < siesta.fdf > siesta.out"
for i, xyz in enumerate(xyz_range):
if xyz < 10:
num = "00000{0}".format(xyz)
elif xyz < 100:
num = "0000{0}".format(xyz)
elif xyz < 1000:
num = "000{0}".format(xyz)
elif xyz < 10000:
num = "00{0}".format(xyz)
else:
raise ValueError("xyz too large?? {0}".format(xyz))
path = "calc_" + num
os.chdir(path)
# Run siesta
subprocess.call("export OMP_NUM_THREADS=1", shell=True)
subprocess.call(siesta_exe, shell=True)
# Run pyscf.nao
subprocess.call("export OMP_NUM_THREADS={0}".format(args_par.np), shell=True)
run_tddft_iter()
os.chdir("../")