ABINIT Documentation is insufficient

Total energy, geometry optimization, DFT+U, spin....

Moderator: bguster

Locked
woffermans
Posts: 41
Joined: Fri Jun 01, 2018 8:22 am
Contact:

ABINIT Documentation is insufficient

Post by woffermans » Thu Jan 31, 2019 10:48 am

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.

woffermans
Posts: 41
Joined: Fri Jun 01, 2018 8:22 am
Contact:

Re: ABINIT Documentation is insufficient

Post by woffermans » Thu Jan 31, 2019 11:27 am

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.

woffermans
Posts: 41
Joined: Fri Jun 01, 2018 8:22 am
Contact:

Re: ABINIT Documentation is insufficient

Post by woffermans » Thu Jan 31, 2019 3:14 pm

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 :)

woffermans
Posts: 41
Joined: Fri Jun 01, 2018 8:22 am
Contact:

Re: ABINIT Documentation is insufficient

Post by woffermans » Mon Feb 04, 2019 2:50 pm

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>

woffermans
Posts: 41
Joined: Fri Jun 01, 2018 8:22 am
Contact:

Re: ABINIT Documentation is insufficient

Post by woffermans » Tue Feb 05, 2019 10:45 am

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>

Locked