Need an example SSCHA run with VASP

[chenyubi14]

Hi,

I understand ASE is supposed to work with VASP. However, I am not able to figure out how to run VASP by editting the example files with quantum espresso. Could you upload a VASP example file?

Thanks!

[diegomartinez2]

Hi,
Using the first tutorial as example.
Just change the calculator from EMT to VASP
From:

 # Import the sscha code
 import sscha, sscha.Ensemble, sscha.SchaMinimizer
 import sscha.Relax, sscha.Utilities
 
 # Import the cellconstructor library to manage phonons
 import cellconstructor as CC, cellconstructor.Phonons
 import cellconstructor.Structure, cellconstructor.calculators
 
 # Import the force field of Gold
 import ase, ase.calculators
 from ase.calculators.emt import EMT
 
 # Import numerical and general pourpouse libraries
 import numpy as np, matplotlib.pyplot as plt
 import sys, os

First change the from ase.calculators.emt import EMT to from ase.calculators.vasp import Vasp.
Later you change calculator = EMT() with calculator = Vasp(prec='Accurate', xc='PBE', ibrion=2, lreal=False) ; now the relax = CC.calculators.Relax(gold_structure, calculator) will use VASP.

Of course this is just an example, you can read more about calling VASP from ASE from here and here.

[chenyubi14]

Thanks a lot! I actually tried to follow calling VASP from ASE, but the ASE manual is very confusing for me.
I ended up trying to use 'Vasp calculator' to read_incar, read_potcar, and read_kpoints. It seems read_potcar('POTCAR') is not a successful setup. I will try your method to directly give the incar parameters. I will let you know about whether it works. Thank you again!

[chenyubi14]

Hi,

I tried exactly the VASP setup you recommended, and I have encountered an error message.

This is the code

import sscha, sscha.Ensemble, sscha.SchaMinimizer, sscha.Relax, sscha.Utilities

import cellconstructor as CC, cellconstructor.Phonons
import cellconstructor.Structure, cellconstructor.calculators

import ase, ase.calculators
from ase.calculators.vasp import Vasp

import numpy as np, matplotlib.pyplot as plt
import sys, os

mater_structure = CC.Structure.Structure()
mater_structure.read_generic_file("POSCAR")

command = '~/bin/vasp.631.pod_std'
calculator = Vasp(
          command=command,
          prec='Accurate',
          xc='PBE',
          ibrion=2,
          encut=400,
          ediff=1e-6,
          nelm=200,
          ismear=0,
          sigma=0.01,
          lreal=False,
          kpts=(1,1,1),
          ncore=2,
          lmaxmix=4,
        )
relax = CC.calculators.Relax(mater_structure, calculator)
mater_structure_relaxed = relax.static_relax()

mater_harmonic_dyn = CC.Phonons.compute_phonons_finite_displacements(
    mater_structure_relaxed, 
    calculator, 
    supercell = (4,4,4))

The last line has generated VASP input files in a new folderdisp_0/. I have no problem running vasp directly by ~/bin/vasp.631.pod_std, but it seems ASE is not happy with the calculator. Below is the error message. Do you know whether my setup is correct? What are the possible problems?

Number of symmetry inequivalent displacements: 1
[pod-login1:89806:0:89806] Caught signal 11 (Segmentation fault: address not mapped to object at address 0x7ffc9b572f40)
==== backtrace ====
0 /usr/lib64/libucs.so.0(+0x17970) [0x7f3aa8af7970]
1 /usr/lib64/libucs.so.0(+0x17b22) [0x7f3aa8af7b22]
2 /home/yubi/bin/vasp.631.pod_std() [0x88d740]
3 /home/yubi/bin/vasp.631.pod_std() [0x104001f]
4 /home/yubi/bin/vasp.631.pod_std() [0x10ddd3c]
5 /home/yubi/bin/vasp.631.pod_std() [0x1bc1624]
6 /home/yubi/bin/vasp.631.pod_std() [0x1b99e2d]
7 /home/yubi/bin/vasp.631.pod_std() [0x4083ee]
8 /lib64/libc.so.6(__libc_start_main+0xf5) [0x7f3aacf84555]
9 /home/yubi/bin/vasp.631.pod_std() [0x4082e9]
===================
Traceback (most recent call last):
File "/home/yubi/work/Au_gold/02_sscha/a5_tut1_vasp_scratch/gold_6_vasp_scratch.py", line 82, in
mater_harmonic_dyn = CC.Phonons.compute_phonons_finite_displacements(mater_structure_relaxed, calculator, supercell = (4,4,4))
File "/home/yubi/repo/miniconda3/envs/sscha/lib/python3.10/site-packages/cellconstructor/Phonons.py", line 4461, in compute_phonons_finite_displacements
return compute_phonons_finite_displacements_sym(structure, ase_calculator, epsilon,
File "/home/yubi/repo/miniconda3/envs/sscha/lib/python3.10/site-packages/cellconstructor/Phonons.py", line 4705, in compute_phonons_finite_displacements_sym
fc = Settings.GoParallel(compute_force, list_of_calculations, reduce_op='+')
File "/home/yubi/repo/miniconda3/envs/sscha/lib/python3.10/site-packages/cellconstructor/Settings.py", line 242, in GoParallel
result += function(x, **kwargs)
File "/home/yubi/repo/miniconda3/envs/sscha/lib/python3.10/site-packages/cellconstructor/Phonons.py", line 4693, in compute_force
energy, forces = calculators.get_energy_forces(ase_calculator, s)
File "/home/yubi/repo/miniconda3/envs/sscha/lib/python3.10/site-packages/cellconstructor/calculators.py", line 56, in get_energy_forces
energy = atm.get_total_energy()
File "/home/yubi/repo/miniconda3/envs/sscha/lib/python3.10/site-packages/ase/atoms.py", line 770, in get_total_energy
return self.get_potential_energy() + self.get_kinetic_energy()
File "/home/yubi/repo/miniconda3/envs/sscha/lib/python3.10/site-packages/ase/atoms.py", line 731, in get_potential_energy
energy = self._calc.get_potential_energy(self)
File "/home/yubi/repo/miniconda3/envs/sscha/lib/python3.10/site-packages/ase/calculators/calculator.py", line 709, in get_potential_energy
energy = self.get_property('energy', atoms)
File "/home/yubi/repo/miniconda3/envs/sscha/lib/python3.10/site-packages/ase/calculators/calculator.py", line 737, in get_property
self.calculate(atoms, [name], system_changes)
File "/home/yubi/repo/miniconda3/envs/sscha/lib/python3.10/site-packages/ase/calculators/vasp/vasp.py", line 325, in calculate
raise calculator.CalculationFailed(
ase.calculators.calculator.CalculationFailed: vasp in disp_0 returned an error: -11

[chenyubi14]

I think I have solve this issue. It is because I'm testing on a login node, which has limited resources. If I change the supercell size to (2,2,2), this problem is gone.
I will keep testing the following tutorial. Thanks!

[chenyubi14]

Hi,
I have encountered another problem regarding VASP.

File "/home/yubi/repo/miniconda3/envs/sscha/lib/python3.10/site-packages/sscha/Cluster.py", line 1572, in
calculators = [cellconstructor_calc.copy() for i in range(0, len(jobs))]
AttributeError: 'Vasp' object has no attribute 'copy'

Basically, the default code assumes quantum espresso calculator obtained from cellconstructor from cellconstructor.calculators import Espresso. However, the vasp calculator is from ASE from ase.calculators.vasp import Vasp. Therefore, 'Vasp' object has no attribute 'copy'. Do you have a way to fix it?

Thanks

[diegomartinez2]

I'm totally not sure about this... but you can define a new function:

    def copy(self):
        """
        Return an identical instance, without inhering the info of the calculation.
        """
        new_class = Vasp(command=command,
          prec='Accurate',
          xc='PBE',
          ibrion=2,
          encut=400,
          ediff=1e-6,
          nelm=200,
          ismear=0,
          sigma=0.01,
          lreal=False,
          kpts=(1,1,1),
          ncore=2,
          lmaxmix=4,)
        return new_class

and then use the setattr() to set the copy function into the object:
setattr(cellconstructor_calc, "copy", copy)
In general it will be:

class MyClass:
    pass

def new_function(self):
    print("This is a new function!")

obj = MyClass()
setattr(obj, "new_function", new_function)

obj.new_function()  # This will call the new function

but I never used this to add a new function, so I do not known if this will work.

EDIT: If this works you have to add the new function just after the creation/definition of the object, as:

calculator_vasp_with_copy = Vasp(command=command,
          prec='Accurate',
          xc='PBE',
          ibrion=2,
          encut=400,
          ediff=1e-6,
          nelm=200,
          ismear=0,
          sigma=0.01,
          lreal=False,
          kpts=(1,1,1),
          ncore=2,
          lmaxmix=4,)
def New_function():
         new_class = Vasp(command=command,
          prec='Accurate',
          xc='PBE',
          ibrion=2,
          encut=400,
          ediff=1e-6,
          nelm=200,
          ismear=0,
          sigma=0.01,
          lreal=False,
          kpts=(1,1,1),
          ncore=2,
          lmaxmix=4,)
        return new_class
setattr(calculator_vasp_with_copy, "copy", New_function)

[chenyubi14]

Hi,

Thanks for the help. After working around this problem, more problems have come out. My current conclusion is that the automatic submission like

ensemble.compute_ensemble(espresso_config, cluster=cluster_config)

and

relax = sscha.Relax.SSCHA(minimizer, calculator,
        N_configs=10,
        max_pop=3,
        save_ensemble=True,
        cluster=cluster_config)

both do not work for VASP.
These functions assumed espresso input format, and using a VASP calculator is not able to write input file.

Thanks. I think I am gonna use the manual submission.

[mesonepigreco]

Hi, unfortunately yes. The automatic submission only works with quantum espresso.
I believe that with some hacks it should be possible to make it run with VASP and others, however the effort may not be worth since you can anyway use the manual submission.
Sooner or later I would like to create a more abstract Cluster class that can work with any calculator. We are also implementing a new Aiida interface, but unfortunately also that one still is program specific as it does not provide a higher abstraction level.
This is on the TODO list for the code.

[diegomartinez2]

Outside of the automatic submission there is this that I use for VASP calculations with Machine Learning. It is a fork of the work of Luigi Ranalli from the University of Vienna.

[chenyubi14]

Thanks.
The package vasp-phonopy-sscha looks useful, but it seems to only work for initial condition.
Please let me know as soon as the Cluster is ready for VASP by commeting on this github issue. I really appreciate it!

[kw600]

I am currently using SSCHA + VASP now. For the mannual submission, one only needs to transform the generated QE inputs into vasp inputs and then run the VASP calculations and finally read the force and energy from OUTCAR. This can be easily done without changing the SSCHA source code. For the automatic submission, one needs to modify the Cluster.py in sscha. The idea is the same as the mannual submission case. The Cluster.py generate submission file for your clusters and one can add cmd lines in the submission file to change the QE input into VASP input and transform the VASP output into QE output.

[Guodonglin-cqust]

I am currently using SSCHA + VASP now. For the mannual submission, one only needs to transform the generated QE inputs into vasp inputs and then run the VASP calculations and finally read the force and energy from OUTCAR. This can be easily done without changing the SSCHA source code. For the automatic submission, one needs to modify the Cluster.py in sscha. The idea is the same as the mannual submission case. The Cluster.py generate submission file for your clusters and one can add cmd lines in the submission file to change the QE input into VASP input and transform the VASP output into QE output.

The atomic positions generated by QE are not in the same order as those in the POSCAR, leading to incorrect forces for atomic positions generated by QE and an increase in free energy rather than a decrease. Do you have any methods to resolve this issue?

[diegomartinez2]

For a work with SrTiO3 (using SSCHA+VASP) I made something that works for that structure but it is not generic.

         for i in range(0, N):
            for j in range(len(type_atoms)):
                for jj in range(0,int(atom_occurrencies[j])):
                    print("[",atom_occurrencies[j]*i+j*N+jj,"]")
                    force =  forces.iloc[atom_occurrencies[j]*i+j*N+jj]
                    force = force.to_frame()
                    force = force.transpose()
                    force.to_csv(f, index=False, header=False, float_format="%16.12f", sep='\t', mode="w+")

Again, this is not a generic method but a patch to reorder the output to VASP. The basic idea is to make a small code that transforms to the POSCAR format of one kind of atom per block.

[Guodonglin-cqust]

For a work with SrTiO3 (using SSCHA+VASP) I made something that works for that structure but it is not generic.

         for i in range(0, N):
            for j in range(len(type_atoms)):
                for jj in range(0,int(atom_occurrencies[j])):
                    print("[",atom_occurrencies[j]*i+j*N+jj,"]")
                    force =  forces.iloc[atom_occurrencies[j]*i+j*N+jj]
                    force = force.to_frame()
                    force = force.transpose()
                    force.to_csv(f, index=False, header=False, float_format="%16.12f", sep='\t', mode="w+")

Again, this is not a generic method but a patch to reorder the output to VASP. The basic idea is to make a small code that transforms to the POSCAR format of one kind of atom per block.

  I have attempted to ensure consistency in atomic order between the POSCAR file and the scf_population.dat file by directly converting the atomic order from POSCAR to scf_population.dat. My initial assumption was that this approach may resolve the discrepancies in atomic ordering.

  However, after executing the python_minimizer.py script, I observed that the free energy continues to increase, indicating that my method of converting the atomic order is not effective. 

[diegomartinez2]

As commented before you can use Luigi Ranalli script to transform from SSCHA to VASP and vice-versa, but probably you need to edit the script to your needs.