Phonon Calculations in CP2K

[Ashley Dickson]

Hello everyone, 

I'm unsure if this is the right place to post this issue as I am using phonopy to generate phonons with CP2K, however I thought I might have a bit more luck finding someone with experience using phonopy specifically with CP2K. The issue I'm having is large negative phonon modes, and large drift force reported by phonopy in the y direction (~0.25 eV/A). I know that I shouldn't be getting negative modes as my structure is certainly dynamically stable. my workflow is as follows:

phonopy --cp2k -d --dim="6 6 2" -c energy.inp 

for i in energy-supercell-*; do
    bn=${i%.*}  # Remove file extension
    num=${bn#*supercell-}  # Extract the number after 'energy-supercell-'

    if [[ -z ${num} ]]; then
        continue
    fi
    dir="$num"

    if [[ -d $dir ]]; then
        echo "Skipping: Directory $dir already exists."
        continue  # Skip to the next iteration instead of exiting
    fi

    mkdir "$dir"

    mv "$i" "$dir/energy.inp"
  #  if [[ "$dir" != "001" ]]; then # Keep Wavefunction file from first simulation to make others quicker
 #       cp 001/MIN-supercell-001-RESTART.wfn "$dir/MIN-supercell-$num-RESTART.wfn"
#    fi

    cd "$dir"
    srun cp2k.psmp -i energy.inp -o result.out
    cd ../

done

phonopy --cp2k -f $(for i in {1..13}; do printf "%03d/MIN-supercell-%03d-forces-1_0.xyz " $i $i; done)

phonopy --cp2k -c energy.inp -p -s --dim="6 6 2" --band " 0.0 0.0 0.0  0.0 0.0 0.5  0.0 0.5 0.5  0.0 0.5 0.0  0.5 0.5 0.0  0.5 0.0 0.0  0.5 0.0 0.5  0.5 0.5 0.5" --mesh 41 41 41

I have tried various supercell sizes, and have found that 4x4x2 should be more than converged enough. My energy.inp file has the coordinates of the relaxed unitcell printed as follows:
    &CELL
      A 3.870961 0 0
      B 0 3.9540005 0
      C 0 0 11.8602555
                PERIODIC XYZ
    &END CELL
   &COORD
Ba 1.9354787033 1.9770037827 9.6771343800
Ba 1.9354787999 1.9770022809 2.1442401159
Y 1.9354802394 1.9770060545 5.8953066752
Cu -0.0000006631 0.0000091908 7.6100822750
Cu -0.0000006426 0.0000055768 4.1724491018
Cu -0.0000078256 0.0000008848 -0.0048129326
O -0.0000004062 1.9770007276 -0.0014139445
O 1.9354801904 0.0000011797 7.3249015553
O 1.9354801707 0.0000009008 4.4580654293
O -0.0000001259 1.9770034877 7.3279693035
O -0.0000001304 1.9770024288 4.4553471390
O -0.0000006180 0.0000004257 9.9643054004
O -0.0000006737 0.0000003017 1.8887411088
      &END COORD

My geometry has been relaxed in a 6x6x2 supercell, with an RMS_FORCE setting of 1e-8 with the CELL_OPT method. The sum of atomic forces at the end lines up with this. I have converged my plane wave cutoff and relative plane wave cutoff (how converged does this need to be?). I worry it may be an issue with negative coordinates in the unit cell, however the periodic XYZ command should wrap these as I understand? 

If anyone has any suggestions as to why my phonons look quite bad I'd be very thankful.

Thanks in advance 

Ashley 

[Krack, Matthias]

Hi Ashley

 

Regarding the CP2K side, you should check if all the finite difference force calculations converged properly to equivalent (magnetic?) states during an atomic displacement. You do not report the applied computational method and I can only guess, but I would have a look at the Cu 3d occupations of the YBCO before and after an atomic displacement.

 

Best

 

Matthias

 

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/cp2k/d5c45f96-776c-4803-90b0-2584e1cd8f4an%40googlegroups.com.

[Ashley Dickson]

Hi Matthias,

thanks for the response. For your information I am using the finite difference method here. Supposing this was the issue, how would you recommend solving it? Presumably I can't explicitly freeze the orbital occupations for the displaced configurations? Perhaps including the wavefunction of the initial geometry as an initial guess might keep the correct spin state.

Also, I wonder if the anisotropic nature of YBCO is the issue here. It is not strictly magnetic, although the copper chains/planes may act as such and mess up the forces. Would it be beneficial to include spin in my calculations maybe? 

Many thanks 

Ashley 

[Ashley Dickson]

Also adding on to my previous comments, I use the following mixing approach:
   &MIXING  T
        METHOD BROYDEN_MIXING
        ALPHA 0.2
        BETA 1.5
        NBUFFER 15
      &END MIXING

Could this be responsible for any differences in convergence issues perhaps?

[Krack, Matthias]

Hi Ashley

 

Even if you could apply easily (which I doubt) a constraint like the freezing of orbital occupations, it would presumable make things worse. I suggest to use the converged wavefunction from the central point as an initial wavefunction for the run with a displaced atom in case you don’t do that already.

 

I guess that the Cu atoms in your system have different oxidation states. Can you identify/assign the oxidation state for each Cu atom before and after displacement?

 

I assumed that you perform UKS runs. If not so far, I would certainly give it a try, because the Cu atoms could exhibit different on-site 3d open-shell configurations while the total spin of the system is still zero. Though that can provide interesting insight, it might also cause new technical (convergence) problems.

 

Best

 

Matthias

 

To view this discussion visit https://groups.google.com/d/msgid/cp2k/d8d75fed-c55e-4288-adff-cdb9e9e31e16n%40googlegroups.com.

[Ashley Dickson]

Hi Matthias,

Thanks you very much for all the help here. In the first instance I will try and start from a converged wavefunction of the perfect cell, hopefully this irons out the issues. I will make sure I print the relevant charge information on my next run to identify if the oxidation states may be playing a role as you mentioned. 

Ideally i'd like to avoid UKS runs primarily because I am validating some results that don't utilise a spin calculation. Although I do believe this may solve my issue, so I may give it a try just to see.

Many thanks for all your help,

Ashley