GAP Usage (Python/Bash)

One needs the following files in the working folder:

INCAR-md, INCAR-static, POSCAR, POTCAR, KPOINTS

#lipai@mail.ustc.edu.cn
import os
MPI="~/bin/mpirun"
VASP="vasp_std"
runvasp=MPI+" -n 64 "+VASP+";"
#runvasp="mpirun -n 24 vasp.5.3.5;"

os.system("mkdir MD;\
           cp POSCAR KPOINTS POTCAR MD;\
           cp INCAR-md MD/INCAR;\
           cd MD;"\
           +runvasp+\
          "outcar2xyz.sh;\
           cp all.xyz ..;\
           cd ..")

gap_fit="gap_fit at_file=all.xyz \
        force_parameter_name=forces \
        gap={ soap cutoff=5.0 \
          atom_sigma=0.5 \
          l_max=6 \
          n_max=10 \
          delta=0.5 \
          covariance_type=dot_product \
          zeta=4 \
          n_sparse=400 \
          add_species=T \
          sparse_method=cur_points} \
e0={C:0:Cu:0} \
default_sigma={0.010 0.20 0.0 0.0} \
do_copy_at_file=F sparse_separate_file=F \
gp_file=gap_soap.xml 2>&1 | grep -v FoX "

os.system(gap_fit)
os.system("cp gap_soap.xml gap_3000.xml")

from ase import Atoms,units
import ase.io
from ase.io.trajectory import Trajectory
from ase.md.nvtberendsen import NVTBerendsen
from quippy.potential import Potential

T=[2000,1800,1600,1400,1200,1000,800,600]
N=[50,  200,  400, 500, 400, 250, 100, 50]

for tt,i in enumerate(T):
    print("GAP MD at "+str(i)+" K")
    calc = Potential(param_filename='gap_soap.xml')
    CuC = ase.io.read("POSCAR")
    CuC.set_calculator(calc)
    dyn = NVTBerendsen(CuC,2*units.fs,temperature=i,taut=0.2*1000*units.fs,fixcm=False)
    traj = Trajectory('gap.traj', 'w', CuC)
    dyn.attach(traj.write, interval=1)
    dyn.run(100*N[tt])
    out_traj = ase.io.read('gap.traj', ':')
    os.system("mkdir "+str(i))
    n=0
    for at in out_traj:
        at.wrap()
        if 'momenta' in at.arrays: del at.arrays['momenta']
        n=n+1
        if n%100==0:
            dir_i=str(i)+"/"+str(n)+"/"
            os.system("mkdir "+dir_i+";\
                       cp KPOINTS POTCAR "+dir_i+";\
                       cp INCAR-static "+dir_i+"INCAR")
            ase.io.write(dir_i+"POSCAR",out_traj[n-1])
            os.system("cd "+dir_i+";"+runvasp+"cd ../..")
    ase.io.write(str(i)+"_gap.xyz",out_traj[99::100])
    os.system("cd "+str(i)+"; \
               for i in *; do \
                 cd $i;\
                 outcar2xyz.sh;\
                 cat all.xyz >> ../../"+str(i)+"_vasp.xyz;\
                 cat all.xyz >> ../../all.xyz;\
                 cd ..;\
               done;\
               cd ..")

    os.system(gap_fit)
    os.system("cp gap_soap.xml gap_"+str(i)+".xml")

To get GP error, use:

quip verbosity=ANALYSIS atoms_filename=test.xyz param_filename=gap_soap.xml E=T F=T L=T calc_args="local_gap_variance" output_file=a.out