ABINIT Discussion Forums

Dear Abinit,

I am trying to perform a berry-phase polarization calculation using the PAW method. Previously I have been able to do this successfully using the FHI pseudopotentials however I would like to know the difference, if any, with a fixed cell/atomic positions, between the two methods. I am currently using the Abinit version 6.0.3 but I have tried in earlier versions also, with no success. After running to SCF convergence I get the message pased below from my log file. Below the log file I supply the input file.

My question is, can one use the berry phase method for polarization calculation with PAW or has this not yet been implemented? Sorry but I have not been able to find it stated explicitly one way or the other in the tutorials. As a follow-up, if it is implemented, does someone know which flags I might have improperly set to perform such a calculation?

Thanks very much for any help with this issue.

Chris

Excerpt from log file:

At SCF step 15, etot is converged :
for the second time, diff in etot= 7.656E-09 < toldfe= 1.000E-07
berryphase_new : enter
berryphase_new : enter
berryphase_new : enter
berryphase_new : enter
berryphase_new : enter
berryphase_new : enter
berryphase_new : enter
berryphase_new : enter
berryphase_new : enter
berryphase_new : enter
berryphase_new : enter
berryphase_new : enter
berryphase_new : enter
berryphase_new : enter
berryphase_new : enter
berryphase_new : enter
Error in cprj_copy: n2 wrong sizes !
Error in cprj_copy: n2 wrong sizes !
Error in cprj_copy: n2 wrong sizes !
MPI: MPI_COMM_WORLD rank 8 has terminated without calling MPI_Finalize()
MPI: aborting job

Input file:

dilatmx 1.1
ecutsm 0.5

#Ground state calculation
toldfe 1.0d-7 # SCF stopping criterion
# causes SCF cycle to stop ions to be moved).
getwfk -1 # This is to speed up the calculation, by restarting
# from previous wavefunctions, transferred from the old
# to the new k-points.
occopt 1
#========================================================
# Printing options
#========================================================
prtwf=1
#========================================================
# Parallelization
#========================================================
#===========================================================
#Setup of SCF calc.
#===========================================================
#Gives the number of band, explicitely (do not take the default)
nband 160 # For an insulator (if described correctly as an insulator
# by DFT), there is no need to include conduction bands
# in response-function calculations

#Exchange-correlation functional
ixc 11 # GGA PBE parameterization
#Definition of the planewave basis set
ecut 20.0 # Maximal kinetic energy cut-off, in Hartree
pawecutdg 40.0
pawovlp -1.0
pawmixdg 1
#===========================================================
#Gives the number of band, explicitely (do not take the default)
#Definition of the k-point grid
#kptopt 2 # For RF Calc.
#kptopt 1 # Automatic generation of k points, taking
kptopt 3 # Automatic generation of k points, no symmetry (eg paw)
kptrlatt 4 0 0 0 4 0 0 0 4
nshiftk 1
shiftk 0.500 0.500 0.500
# shiftk 0.000 0.000 0.000
#===========================================================
#Definition of the SCF procedure
iprcel 49 # compute RPA dielectric matrix at first step, recompute at later step, take into account mixing factor diemix
diemix 0.50 # factor of the preconditioned residual potential to be transferred in the SCF cycle
iscf 17 # Anderson mixing of potential/Pullay mixing with prev. 2 iter
nstep 80 # Maximal number of SCF cycles
diemac 3.0 # Although this is not mandatory, it is worth to
# precondition the SCF cycle. The model dielectric

#===========================================================
# Berry phase calculation of the polarization
#===========================================================
berryopt -1
rfdir 1 1 1
nbdbuf 0
#===========================================================
#######################################################################
#Setup of the cell

#Definition of the unit cell
I leave the cell structure out, but I can send it if it's helpful

This has certainly been coded for PAW, although I am not certain what is the status in 6.0.3. Try running a serial calculation (eventually with low cutoff for testing). Otherwise your input looks ok

Error in cprj_copy: n2 wrong sizes !

smacks of a bug. It looks like the 6.1.2 code is more coherent, and there may have been a bug in 6.0.3:

Check line 273 in file src/79_seqpar_mpi/berryphase_new.F90 and print out dtefield%fnkpt, nkpt, and mkmem. If the last 2 are different that's the bug.

1) Did you use mkmem 0?
2) check that your k-grid shift is consistent with the full symmetry of the cell (basically must be cubic for 0.5 0.5 0.5). For hexagonal or tetragonal, use 0 0 0.5 and otherwise try with 0 0 0. We at least need the unit cell structure, if not the positions.

If we can confirm that there is a real bug, we can try to fix it in 6.0.4 (MG?)

Matthieu

I will post your response in quotes below. Mine are not in quotes.

"Check line 273 in file src/79_seqpar_mpi/berryphase_new.F90 and print out dtefield%fnkpt, nkpt, and mkmem. If the last 2 are different that's the bug."

Inside the log file I find the line:

invars1: mkmem undefined in the input file. Use default mkmem = nkpt
invars1: With nkpt_me= 4 and mkmem = 64, ground state wf handled in core.

I do not know if this is the same information as you recommend in the above quote. I am trying a different shiftk at the moment. If that doesn't work I will go back and print out the above information.

"1) Did you use mkmem 0?"

I did not set mkmem, but according to the Abinit web page input variables section the default is mkmem=nkpt.

"2) check that your k-grid shift is consistent with the full symmetry of the cell (basically must be cubic for 0.5 0.5 0.5). For hexagonal or tetragonal, use 0 0 0.5 and otherwise try with 0 0 0. We at least need the unit cell structure, if not the positions."

I have a tetragonal cell:
#Definition of the unit cell
acell 1.4603809311E+01 1.4603809311E+01 1.5250022218E+01
rprim= 1.0000000000E+00 0.0000000000E+00 0.0000000000E+00
0.0000000000E+00 1.0000000000E+00 0.0000000000E+00
0.0000000000E+00 0.0000000000E+00 1.0000000000E+00

Sorry for leaving that information out in my previous e-mail and thanks for recognizing it. I am now hoping that my problem is an improper choice of shiftk. I used the 0.5 0.5 0.5. I am now trying the shiftk = 0.0 0.0 0.5. I will keep you posted.

Thanks very much for your reply.

Hi, I am the one currently doing the electric field + PAW developments. In 6.0 the serial calculation works well but not the MPI version. The 6.1 branch is currently developing and improving rapidly. There will be major updates to electric field plus PAW shortly, probably within 2 weeks.

Joe

I have tried using both shiftk=0 0 0 and 0 0 .5. Both crash with the same error message. I am also running it in serial mode but this will take a few days. I will post the result as to whether it ran successfully or not.

I have tried a smaller system (rhombohedral BTO with a 5 atom cell). The serial version works fine and I get a reasonable output. The parallel version crashes at about the same point (shortly after entering berryphase_new--although I don't know if it's still in there when it crashes).

I would gladly post my input file and error message for this if it is of interest.

Thanks again to both of you for your comments.

Just want to finish this posting by stating that it worked in serial. I am able to do the parallel berry phase calculation using the abinit 6.2.2.

Thanks for the wonderful, well managed code.