Dear ABINIT friends,
I'm currently performing STRING calculations.
There is documentation to this topic within the ABINIT documentation set:
https://docs.abinit.org/topics/TransPath/
However, after reading the above, I'm still left behind with questions:
1) How do I restart a STRING calculation?
2) How do I ensure that the previous WFK files are reused in the new run?
At the moment there are some bugs in ABINIT, which prevent a swift calculation. A string calculation can in principle be restarted, but the calculation stops just before finalising the TIME STEP. A segmentation error occurs.
So I want to extract the atom positions of the images out of the HIST.nc file. There is no documentation on how to do this. At least I cannot find it.
3) How can atom positions of the images being extracted from HIST.nc file?
I would love to contribute to the ABINIT documentation, but I will firstly need answers to my questions.
So if someone is willing to respond, I would be very grateful.
ABINIT Documentation is insufficient
Moderator: bguster
-
- Posts: 41
- Joined: Fri Jun 01, 2018 8:22 am
- Contact:
-
- Posts: 41
- Joined: Fri Jun 01, 2018 8:22 am
- Contact:
Re: ABINIT Documentation is insufficient
Dear ABINIT friends,
To answer my questions partially, it is possible to dump the content of the NetCDF file (HIST.nc) into ASCII.
ncdump is your friend:
https://www.unidata.ucar.edu/software/n ... II-or-text
I don't know how to get ncdump. It was installed on my macOS:
/Users/willem/anaconda3/bin/ncdump
So it came with some package within anaconda.
At least one can now read the content of the HIST.nc file.
This is exactly the documentation, that I'm missing.
To answer my questions partially, it is possible to dump the content of the NetCDF file (HIST.nc) into ASCII.
ncdump is your friend:
https://www.unidata.ucar.edu/software/n ... II-or-text
I don't know how to get ncdump. It was installed on my macOS:
/Users/willem/anaconda3/bin/ncdump
So it came with some package within anaconda.
At least one can now read the content of the HIST.nc file.
This is exactly the documentation, that I'm missing.
-
- Posts: 41
- Joined: Fri Jun 01, 2018 8:22 am
- Contact:
Re: ABINIT Documentation is insufficient
Dear ABINIT friends,
It is relatively easy to extract the atom positions of the images of the STRING calculation, stored in the HIST.nc file.
I hope you know a bit of bash scripting:
<snip>
#! /bin/sh
###############################################################################
#
#
# a quick hack to extract Cartesian atom positions from a
# STRING calculation, stored in a NetCDF history file
#
#
# Willem.Offermans@vito.be 31/01/2019
#
# Dump content of HIST.nc file into a ASCII text file:
# for example:
# `ncdump pdo_HIST.nc > pdo_HIST.txt`
#
#
###############################################################################
HISTFILE='./pdo_HIST.txt'
NoAtomsPerImage=`grep "natom =" ${HISTFILE} | awk '{print $3}'`
NoImages=`grep "nimage =" ${HISTFILE} | awk '{print $3}'`
xcartPOS=`grep -n 'xcart =' ${HISTFILE} | awk 'BEGIN{FS=":"}{print $1}'`
xredPOS=`grep -n 'xred =' ${HISTFILE} | awk 'BEGIN{FS=":"}{print $1}'`
xredLinesPos=$((${xredPOS}-${xcartPOS}-2))
grep -A${xredLinesPos} "xcart =" ${HISTFILE} | \
tail -n${xredLinesPos} | \
awk '
# set some variables for awk
BEGIN{lineCounter=0;atomCounter=0;imageCounter=1;tstep=1}
{
# reset image lineCounter to 1, if all images of a TIME STEP have been processed
if((imageCounter-1) % '${NoImages}' ==0)
{
imageCounter=1;
};
# print TIME STEP if lineCounter is multiple of timeStep*NumberOfImages*NumerOfAtomsPerImage
# update timeStep
if((lineCounter % (tstep * '${NoImages}' * '${NoAtomsPerImage}')) == 0)
{
if(lineCounter!=0){tstep++};
print "================================================================================";
printf "STRING METHOD - TIME STEP %.i\n",tstep;
print "================================================================================";
};
# print xangst for each image
# reset atomCounter
# update imageCounter
if((atomCounter % '${NoAtomsPerImage}')==0)
{
printf "xangst_%iimg\n", imageCounter;
atomCounter=0;
imageCounter++
};
{printf"%.9f %.9f %.9f \n", $1*0.52917721067,$2*0.52917721067,$3*0.52917721067};
atomCounter++;
lineCounter++;
} '
exit 0
</snip>
I hope it will help you as well
It is relatively easy to extract the atom positions of the images of the STRING calculation, stored in the HIST.nc file.
I hope you know a bit of bash scripting:
<snip>
#! /bin/sh
###############################################################################
#
#
# a quick hack to extract Cartesian atom positions from a
# STRING calculation, stored in a NetCDF history file
#
#
# Willem.Offermans@vito.be 31/01/2019
#
# Dump content of HIST.nc file into a ASCII text file:
# for example:
# `ncdump pdo_HIST.nc > pdo_HIST.txt`
#
#
###############################################################################
HISTFILE='./pdo_HIST.txt'
NoAtomsPerImage=`grep "natom =" ${HISTFILE} | awk '{print $3}'`
NoImages=`grep "nimage =" ${HISTFILE} | awk '{print $3}'`
xcartPOS=`grep -n 'xcart =' ${HISTFILE} | awk 'BEGIN{FS=":"}{print $1}'`
xredPOS=`grep -n 'xred =' ${HISTFILE} | awk 'BEGIN{FS=":"}{print $1}'`
xredLinesPos=$((${xredPOS}-${xcartPOS}-2))
grep -A${xredLinesPos} "xcart =" ${HISTFILE} | \
tail -n${xredLinesPos} | \
awk '
# set some variables for awk
BEGIN{lineCounter=0;atomCounter=0;imageCounter=1;tstep=1}
{
# reset image lineCounter to 1, if all images of a TIME STEP have been processed
if((imageCounter-1) % '${NoImages}' ==0)
{
imageCounter=1;
};
# print TIME STEP if lineCounter is multiple of timeStep*NumberOfImages*NumerOfAtomsPerImage
# update timeStep
if((lineCounter % (tstep * '${NoImages}' * '${NoAtomsPerImage}')) == 0)
{
if(lineCounter!=0){tstep++};
print "================================================================================";
printf "STRING METHOD - TIME STEP %.i\n",tstep;
print "================================================================================";
};
# print xangst for each image
# reset atomCounter
# update imageCounter
if((atomCounter % '${NoAtomsPerImage}')==0)
{
printf "xangst_%iimg\n", imageCounter;
atomCounter=0;
imageCounter++
};
{printf"%.9f %.9f %.9f \n", $1*0.52917721067,$2*0.52917721067,$3*0.52917721067};
atomCounter++;
lineCounter++;
} '
exit 0
</snip>
I hope it will help you as well

-
- Posts: 41
- Joined: Fri Jun 01, 2018 8:22 am
- Contact:
Re: ABINIT Documentation is insufficient
Dear ABINIT friends,
For the people who like to work with python. The same script as above, but now in python code:
<snip>
#! /Users/willem/anaconda3/bin/python3.7
###############################################################################
#
#
# a quick hack to extract Cartesian atom positions from a
# STRING calculation, stored in a NetCDF history file
#
#
# Willem.Offermans@vito.be 04/02/2019
#
#
#
###############################################################################
from netCDF4 import Dataset
import numpy as np
# Read dataset from NetCDF file
myDS = Dataset("pdo_HIST.nc", "r")
# recalculate to Angstrom
timeStepsImagesAtomCartPos=myDS.variables['xcart'][:].data*0.5291772105638411
for timeStepCounter in range(0,len(timeStepsImagesAtomCartPos)):
print('================================================================================')
print('STRING METHOD - TIME STEP '+str(timeStepCounter+1))
print('================================================================================')
for imageCounter in range(0,len(timeStepsImagesAtomCartPos[timeStepCounter])):
print('xangst_'+str(imageCounter+1)+'img')
print('\n'.join(("%20.14f "*timeStepsImagesAtomCartPos[timeStepCounter][imageCounter].shape[1] % tuple(x) for x in timeStepsImagesAtomCartPos[timeStepCounter][imageCounter])))
</snip>
For the people who like to work with python. The same script as above, but now in python code:
<snip>
#! /Users/willem/anaconda3/bin/python3.7
###############################################################################
#
#
# a quick hack to extract Cartesian atom positions from a
# STRING calculation, stored in a NetCDF history file
#
#
# Willem.Offermans@vito.be 04/02/2019
#
#
#
###############################################################################
from netCDF4 import Dataset
import numpy as np
# Read dataset from NetCDF file
myDS = Dataset("pdo_HIST.nc", "r")
# recalculate to Angstrom
timeStepsImagesAtomCartPos=myDS.variables['xcart'][:].data*0.5291772105638411
for timeStepCounter in range(0,len(timeStepsImagesAtomCartPos)):
print('================================================================================')
print('STRING METHOD - TIME STEP '+str(timeStepCounter+1))
print('================================================================================')
for imageCounter in range(0,len(timeStepsImagesAtomCartPos[timeStepCounter])):
print('xangst_'+str(imageCounter+1)+'img')
print('\n'.join(("%20.14f "*timeStepsImagesAtomCartPos[timeStepCounter][imageCounter].shape[1] % tuple(x) for x in timeStepsImagesAtomCartPos[timeStepCounter][imageCounter])))
</snip>
-
- Posts: 41
- Joined: Fri Jun 01, 2018 8:22 am
- Contact:
Re: ABINIT Documentation is insufficient
Dear ABINIT friends,
For the people who like to work with python and ASE, the following script will convert an ABINIT NetCDF History file into an ASE trajectory file:
<snip>
#! /Users/willem/anaconda3/bin/python3.7
###############################################################################
#
#
# Create a ASE trajectory file from ABINIT NetCDF history file
#
#
# Willem.Offermans@vito.be 04/02/2019
#
#
#
###############################################################################
# Don't reinvent the wheel
# but import the good stuff others have made
import datetime
import numpy as np
import os
import shutil
import sys
from ase.calculators.singlepoint import SinglePointCalculator
from ase import Atoms
from ase.io.trajectory import Trajectory
from ase.units import create_units
from netCDF4 import Dataset
# The program needs an ABINIT HISTORY NetCDF file as input
if len(sys.argv) != 2:
print('usage: ' + str(sys.argv[0]) + ' <ABINIT HISTORY NetCDF FILE>')
exit(0)
# Set an output file
ABINITtrajectoryFile = './ABINIT.traj'
# CODATA 2014
units = create_units('2014')
# Read dataset from NetCDF file/ABINIT HISTORY file
myDS = Dataset(str(sys.argv[1]), "r")
# List of integers of considered nuclei
znuclList = list(map(int,myDS.variables['znucl'][:].data.tolist()))
# Determine the length of znuclList
znuclLength = len(znuclList)
# Determine the length of typat
typatLength = len(myDS.variables['typat'][:])
# Create a typatLength x znuclLength matrix with only zeros
znuclMask = np.zeros( (typatLength, znuclLength) )
# The matrix will serve as a mask for the nuclei, once it is
# filled with ones at the right positions
# Set the column with index i to 1 by traversing the rows of
# the znuclMask. The index i is from the typat vector
for typAtomCounter in range(0,len(myDS.variables['typat'][:].data)):
columnIndex = int(myDS.variables['typat'][typAtomCounter].data.tolist())-1
znuclMask[typAtomCounter,columnIndex]=1
timeStepsImagesAtomCartPos=myDS.variables['xcart'][:].data*units.Bohr
# Ensure that ABINIT trajectory file is not overwritten
suffix = datetime.datetime.now().strftime("%y%m%d_%H%M%S")
if os.path.exists(ABINITtrajectoryFile):
shutil.copy2(ABINITtrajectoryFile, ABINITtrajectoryFile + '_' + suffix)
os.remove(ABINITtrajectoryFile)
# Open ABINIT trajectory file to write data
nebTraj = Trajectory(ABINITtrajectoryFile, 'a')
# Loop over TIME STEPS
for timeStepCounter in range(0,len(timeStepsImagesAtomCartPos)):
# Loop over images
for imageCounter in range(0,len(timeStepsImagesAtomCartPos[timeStepCounter])):
image=Atoms(
symbols=list(map(int,np.matmul(znuclMask,znuclList))),
positions=timeStepsImagesAtomCartPos[timeStepCounter][imageCounter],
cell=myDS.variables['rprimd'][timeStepCounter][imageCounter].data*units.Bohr,
pbc=[1, 1, 0]
)
image.set_calculator(SinglePointCalculator(
atoms=image,
energy=myDS.variables['etotal'][timeStepCounter][imageCounter]*units.Hartree,
forces=myDS.variables['fcart'][timeStepCounter][imageCounter].data*units.Hartree/units.Bohr)
)
nebTraj.write(image)
nebTraj.close()
</snip>
For the people who like to work with python and ASE, the following script will convert an ABINIT NetCDF History file into an ASE trajectory file:
<snip>
#! /Users/willem/anaconda3/bin/python3.7
###############################################################################
#
#
# Create a ASE trajectory file from ABINIT NetCDF history file
#
#
# Willem.Offermans@vito.be 04/02/2019
#
#
#
###############################################################################
# Don't reinvent the wheel
# but import the good stuff others have made
import datetime
import numpy as np
import os
import shutil
import sys
from ase.calculators.singlepoint import SinglePointCalculator
from ase import Atoms
from ase.io.trajectory import Trajectory
from ase.units import create_units
from netCDF4 import Dataset
# The program needs an ABINIT HISTORY NetCDF file as input
if len(sys.argv) != 2:
print('usage: ' + str(sys.argv[0]) + ' <ABINIT HISTORY NetCDF FILE>')
exit(0)
# Set an output file
ABINITtrajectoryFile = './ABINIT.traj'
# CODATA 2014
units = create_units('2014')
# Read dataset from NetCDF file/ABINIT HISTORY file
myDS = Dataset(str(sys.argv[1]), "r")
# List of integers of considered nuclei
znuclList = list(map(int,myDS.variables['znucl'][:].data.tolist()))
# Determine the length of znuclList
znuclLength = len(znuclList)
# Determine the length of typat
typatLength = len(myDS.variables['typat'][:])
# Create a typatLength x znuclLength matrix with only zeros
znuclMask = np.zeros( (typatLength, znuclLength) )
# The matrix will serve as a mask for the nuclei, once it is
# filled with ones at the right positions
# Set the column with index i to 1 by traversing the rows of
# the znuclMask. The index i is from the typat vector
for typAtomCounter in range(0,len(myDS.variables['typat'][:].data)):
columnIndex = int(myDS.variables['typat'][typAtomCounter].data.tolist())-1
znuclMask[typAtomCounter,columnIndex]=1
timeStepsImagesAtomCartPos=myDS.variables['xcart'][:].data*units.Bohr
# Ensure that ABINIT trajectory file is not overwritten
suffix = datetime.datetime.now().strftime("%y%m%d_%H%M%S")
if os.path.exists(ABINITtrajectoryFile):
shutil.copy2(ABINITtrajectoryFile, ABINITtrajectoryFile + '_' + suffix)
os.remove(ABINITtrajectoryFile)
# Open ABINIT trajectory file to write data
nebTraj = Trajectory(ABINITtrajectoryFile, 'a')
# Loop over TIME STEPS
for timeStepCounter in range(0,len(timeStepsImagesAtomCartPos)):
# Loop over images
for imageCounter in range(0,len(timeStepsImagesAtomCartPos[timeStepCounter])):
image=Atoms(
symbols=list(map(int,np.matmul(znuclMask,znuclList))),
positions=timeStepsImagesAtomCartPos[timeStepCounter][imageCounter],
cell=myDS.variables['rprimd'][timeStepCounter][imageCounter].data*units.Bohr,
pbc=[1, 1, 0]
)
image.set_calculator(SinglePointCalculator(
atoms=image,
energy=myDS.variables['etotal'][timeStepCounter][imageCounter]*units.Hartree,
forces=myDS.variables['fcart'][timeStepCounter][imageCounter].data*units.Hartree/units.Bohr)
)
nebTraj.write(image)
nebTraj.close()
</snip>