add bcc interstitial calculation and change some default parameters in lmp

dpgen/auto_test/Interstitial.py

@@ -172,6 +172,120 @@ def make_confs(self,
                     dumpfn(self.supercell, 'supercell.json')
                 os.chdir(cwd)
 
+
+                if 'bcc_self' in self.parameter and self.parameter['bcc_self']:
+                    os.chdir(path_to_work)
+                    with open('POSCAR', 'r') as fin:
+                        fin.readline()
+                        scale = float(fin.readline().split()[0])
+                        latt_param = float(fin.readline().split()[0])
+                        latt_param *= scale
+
+                    if not os.path.isfile('task.000000/POSCAR'):
+                        raise RuntimeError("need task.000000 structure as reference")
+
+                    with open('task.000000/POSCAR', 'r') as fin:
+                        pos_line = fin.read().split('\n')
+
+                    super_latt_param = float(pos_line[2].split()[0])
+
+                    output_task1 = os.path.join(path_to_work, 'task.%06d' % (len(dss)))
+                    os.makedirs(output_task1, exist_ok=True)
+                    os.chdir(output_task1)
+                    task_list.append(output_task1)
+                    with open(insert_element_task, 'a+') as fout:
+                        print(self.insert_ele[0], file=fout)
+                    dumpfn(self.supercell, 'supercell.json')
+                    pos_line[-2] = '%.6f' % float(latt_param/4/super_latt_param) + ' ' + '%.6f' % float(latt_param/2/super_latt_param) + ' 0.000000 ' + self.insert_ele[0]
+                    with open('POSCAR', 'w+') as fout:
+                        for ii in pos_line:
+                                print(ii, file=fout)
+                    print('gen bcc tetrahedral')
+                    os.chdir(cwd)
+
+                    output_task2 = os.path.join(path_to_work, 'task.%06d' % (len(dss)+1))
+                    os.makedirs(output_task2, exist_ok=True)
+                    os.chdir(output_task2)
+                    task_list.append(output_task2)
+                    with open(insert_element_task, 'a+') as fout:
+                        print(self.insert_ele[0], file=fout)
+                    dumpfn(self.supercell, 'supercell.json')
+                    pos_line[-2] = '%.6f' % float(latt_param/2/super_latt_param) + ' ' + '%.6f' % float(latt_param/2/super_latt_param) + ' 0.000000 ' + self.insert_ele[0]
+                    with open('POSCAR', 'w+') as fout:
+                        for ii in pos_line:
+                                print(ii, file=fout)
+                    print('gen bcc octahedral')
+                    os.chdir(cwd)
+
+                    output_task3 = os.path.join(path_to_work, 'task.%06d' % (len(dss)+2))
+                    os.makedirs(output_task3, exist_ok=True)
+                    os.chdir(output_task3)
+                    task_list.append(output_task3)
+                    with open(insert_element_task, 'a+') as fout:
+                        print(self.insert_ele[0], file=fout)
+                    dumpfn(self.supercell, 'supercell.json')
+                    pos_line[-2] = '%.6f' % float(latt_param/4/super_latt_param) + ' ' + '%.6f' % float(latt_param/4/super_latt_param) + ' ' + '%.6f' % float(latt_param/4/super_latt_param) + ' ' + self.insert_ele[0]
+                    with open('POSCAR', 'w+') as fout:
+                        for ii in pos_line:
+                                print(ii, file=fout)
+                    print('gen bcc crowdion')
+                    os.chdir(cwd)
+
+                    for idx, ii in enumerate(pos_line):
+                        ss = ii.split()
+                        if len(ss) > 3:
+                            if abs(latt_param/2/super_latt_param - float(ss[0])) < 1e-5 and abs(latt_param/2/super_latt_param - float(ss[1])) < 1e-5 and abs(latt_param/2/super_latt_param - float(ss[2])) < 1e-5:
+                                replace_label = idx
+
+                    output_task4 = os.path.join(path_to_work, 'task.%06d' % (len(dss)+3))
+                    os.makedirs(output_task4, exist_ok=True)
+                    os.chdir(output_task4)
+                    task_list.append(output_task4)
+                    with open(insert_element_task, 'a+') as fout:
+                        print(self.insert_ele[0], file=fout)
+                    dumpfn(self.supercell, 'supercell.json')
+                    pos_line[-2] = '%.6f' % float(latt_param/3/super_latt_param) + ' ' + '%.6f' % float(latt_param/3/super_latt_param) + ' ' + '%.6f' % float(latt_param/3/super_latt_param) + ' ' + self.insert_ele[0]
+                    pos_line[replace_label] = '%.6f' % float(latt_param/3*2/super_latt_param) + ' ' + '%.6f' % float(latt_param/3*2/super_latt_param) + ' ' + '%.6f' % float(latt_param/3*2/super_latt_param) + ' ' + self.insert_ele[0]
+
+                    with open('POSCAR', 'w+') as fout:
+                        for ii in pos_line:
+                                print(ii, file=fout)
+                    print('gen bcc <111> dumbbell')
+                    os.chdir(cwd)
+
+                    output_task5 = os.path.join(path_to_work, 'task.%06d' % (len(dss)+4))
+                    os.makedirs(output_task5, exist_ok=True)
+                    os.chdir(output_task5)
+                    task_list.append(output_task5)
+                    with open(insert_element_task, 'a+') as fout:
+                        print(self.insert_ele[0], file=fout)
+                    dumpfn(self.supercell, 'supercell.json')
+                    pos_line[-2] = '%.6f' % float((latt_param+2.1/2**0.5)/2/super_latt_param) + ' ' + '%.6f' % float((latt_param-2.1/2**0.5)/2/super_latt_param) + ' ' + '%.6f' % float(latt_param/2/super_latt_param) + ' ' + self.insert_ele[0]
+                    pos_line[replace_label] = '%.6f' % float((latt_param-2.1/2**0.5)/2/super_latt_param) + ' ' + '%.6f' % float((latt_param+2.1/2**0.5)/2/super_latt_param) + ' ' + '%.6f' % float(latt_param/2/super_latt_param) + ' ' + self.insert_ele[0]
+
+                    with open('POSCAR', 'w+') as fout:
+                        for ii in pos_line:
+                                print(ii, file=fout)
+                    print('gen bcc <110> dumbbell')
+                    os.chdir(cwd)
+
+                    output_task6 = os.path.join(path_to_work, 'task.%06d' % (len(dss)+5))
+                    os.makedirs(output_task6, exist_ok=True)
+                    os.chdir(output_task6)
+                    task_list.append(output_task6)
+                    with open(insert_element_task, 'a+') as fout:
+                        print(self.insert_ele[0], file=fout)
+                    dumpfn(self.supercell, 'supercell.json')
+                    pos_line[-2] = '%.6f' % float(latt_param/2/super_latt_param) + ' ' + '%.6f' % float(latt_param/2/super_latt_param) + ' ' + '%.6f' % float((latt_param-2.1)/2/super_latt_param) + ' ' + self.insert_ele[0]
+                    pos_line[replace_label] = '%.6f' % float(latt_param/2/super_latt_param) + ' ' + '%.6f' % float(latt_param/2/super_latt_param) + ' ' + '%.6f' % float((latt_param+2.1)/2/super_latt_param) + ' ' + self.insert_ele[0]
+
+                    with open('POSCAR', 'w+') as fout:
+                        for ii in pos_line:
+                                print(ii, file=fout)
+                    print('gen bcc <100> dumbbell')
+                    os.chdir(cwd)
+
+
         return task_list
 
     def post_process(self, task_list):

dpgen/auto_test/common_equi.py

@@ -39,30 +39,34 @@ def make_equi(confs,
     # ...
     cwd = os.getcwd()
     # generate poscar for single element crystal
-    if len(ele_list) == 1:
+    if len(ele_list) == 1 or 'single' in inter_param:
+        if 'single' in inter_param:
+            element_label = int(inter_param['single'])
+        else:
+            element_label = 0
         for ii in conf_dirs:
             os.chdir(ii)
             crys_type = ii.split('/')[-1]
             dlog.debug('crys_type: %s' % crys_type)
             dlog.debug('pwd: %s' % os.getcwd())
             if crys_type == 'std-fcc':
                 if not os.path.exists('POSCAR'):
-                    crys.fcc1(ele_list[0]).to('POSCAR', 'POSCAR')
+                    crys.fcc1(ele_list[element_label]).to('POSCAR', 'POSCAR')
             elif crys_type == 'std-hcp':
                 if not os.path.exists('POSCAR'):
-                    crys.hcp(ele_list[0]).to('POSCAR', 'POSCAR')
+                    crys.hcp(ele_list[element_label]).to('POSCAR', 'POSCAR')
             elif crys_type == 'std-dhcp':
                 if not os.path.exists('POSCAR'):
-                    crys.dhcp(ele_list[0]).to('POSCAR', 'POSCAR')
+                    crys.dhcp(ele_list[element_label]).to('POSCAR', 'POSCAR')
             elif crys_type == 'std-bcc':
                 if not os.path.exists('POSCAR'):
-                    crys.bcc(ele_list[0]).to('POSCAR', 'POSCAR')
+                    crys.bcc(ele_list[element_label]).to('POSCAR', 'POSCAR')
             elif crys_type == 'std-diamond':
                 if not os.path.exists('POSCAR'):
-                    crys.diamond(ele_list[0]).to('POSCAR', 'POSCAR')
+                    crys.diamond(ele_list[element_label]).to('POSCAR', 'POSCAR')
             elif crys_type == 'std-sc':
                 if not os.path.exists('POSCAR'):
-                    crys.sc(ele_list[0]).to('POSCAR', 'POSCAR')
+                    crys.sc(ele_list[element_label]).to('POSCAR', 'POSCAR')
 
             os.chdir(cwd)
     task_dirs = []

dpgen/auto_test/lib/lammps.py

@@ -115,7 +115,7 @@ def inter_deepmd(param):
             ret += "%s out_freq 10 out_file model_devi.out\n" % model_list
         else:
             ret += models[0] + '\n'
-    ret += "pair_coeff\n"
+    ret += "pair_coeff * *\n"
     return ret
 
 
@@ -175,7 +175,7 @@ def make_lammps_eval(conf, type_map, interaction, param):
     ret += "clear\n"
     ret += "units 	metal\n"
     ret += "dimension	3\n"
-    ret += "boundary	p	p    p\n"
+    ret += "boundary	p p p\n"
     ret += "atom_style	atomic\n"
     ret += "box         tilt large\n"
     ret += "read_data   %s\n" % conf
@@ -212,7 +212,7 @@ def make_lammps_eval(conf, type_map, interaction, param):
 
 
 def make_lammps_equi(conf, type_map, interaction, param,
-                     etol=1e-12, ftol=1e-6,
+                     etol=0, ftol=1e-10,
                      maxiter=5000, maxeval=500000,
                      change_box=True):
     type_map_list = element_list(type_map)
@@ -224,7 +224,7 @@ def make_lammps_equi(conf, type_map, interaction, param,
     ret += "clear\n"
     ret += "units 	metal\n"
     ret += "dimension	3\n"
-    ret += "boundary	p	p    p\n"
+    ret += "boundary	p p p\n"
     ret += "atom_style	atomic\n"
     ret += "box         tilt large\n"
     ret += "read_data   %s\n" % conf
@@ -268,7 +268,7 @@ def make_lammps_equi(conf, type_map, interaction, param,
 
 
 def make_lammps_elastic(conf, type_map, interaction, param,
-                        etol=1e-12, ftol=1e-6,
+                        etol=0, ftol=1e-10,
                         maxiter=5000, maxeval=500000):
     type_map_list = element_list(type_map)
 
@@ -279,7 +279,7 @@ def make_lammps_elastic(conf, type_map, interaction, param,
     ret += "clear\n"
     ret += "units 	metal\n"
     ret += "dimension	3\n"
-    ret += "boundary	p	p    p\n"
+    ret += "boundary	p p p\n"
     ret += "atom_style	atomic\n"
     ret += "box         tilt large\n"
     ret += "read_data   %s\n" % conf
@@ -314,7 +314,7 @@ def make_lammps_elastic(conf, type_map, interaction, param,
 
 def make_lammps_press_relax(conf, type_map, scale2equi, interaction, param,
                             B0=70, bp=0,
-                            etol=1e-12, ftol=1e-6,
+                            etol=0, ftol=1e-10,
                             maxiter=5000, maxeval=500000):
     type_map_list = element_list(type_map)
 
@@ -371,7 +371,7 @@ def make_lammps_press_relax(conf, type_map, scale2equi, interaction, param,
 
 
 def make_lammps_phonon(conf, masses, interaction, param,
-                       etol=1e-12, ftol=1e-6,
+                       etol=0, ftol=1e-10,
                        maxiter=5000, maxeval=500000):
     """
     make lammps input for elastic calculation
@@ -380,7 +380,7 @@ def make_lammps_phonon(conf, masses, interaction, param,
     ret += "clear\n"
     ret += "units 	metal\n"
     ret += "dimension	3\n"
-    ret += "boundary	p	p    p\n"
+    ret += "boundary	p p p\n"
     ret += "atom_style	atomic\n"
     ret += "box         tilt large\n"
     ret += "read_data   %s\n" % conf

tests/auto_test/equi/vasp/CONTCAR_V_bcc

@@ -0,0 +1,13 @@
+V2                                      
+   1.0000000000000000     
+     2.9975462819866885    0.0000000000000000    0.0000000000000000
+     0.0000000000000000    2.9975462819866885    0.0000000000000000
+     0.0000000000000000    0.0000000000000000    2.9975462819866885
+   V 
+     2
+Direct
+  0.0000000000000000  0.0000000000000000  0.0000000000000000
+  0.5000000000000000  0.5000000000000000  0.5000000000000000
+
+  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00

tests/auto_test/test_interstitial.py

@@ -0,0 +1,101 @@
+import os, sys, json, glob, shutil
+import dpdata
+import numpy as np
+import unittest
+import dpdata
+from monty.serialization import loadfn, dumpfn
+from pymatgen.core import Structure
+from pymatgen.io.vasp import Incar
+from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
+from pymatgen.analysis.defects.core import Interstitial as pmg_Interstitial
+
+sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))
+__package__ = 'auto_test'
+
+from .context import make_kspacing_kpoints
+from .context import setUpModule
+
+from dpgen.auto_test.Interstitial import Interstitial
+
+
+class TestInterstitial(unittest.TestCase):
+
+    def setUp(self):
+        _jdata = {
+            "structures": ["confs/std-bcc"],
+            "interaction": {
+                "type": "vasp",
+                "incar": "vasp_input/INCAR.rlx",
+                "potcar_prefix": "vasp_input",
+                "potcars": {"V": "POTCAR"}
+            },
+            "properties": [
+                {
+                    "type": "interstitial",
+                    "supercell": [1, 1, 1],
+                    "insert_ele": ["V"],
+                    "bcc_self": True
+                }
+            ]
+        }
+
+        self.equi_path = 'confs/std-bcc/relaxation/relax_task'
+        self.source_path = 'equi/vasp'
+        self.target_path = 'confs/std-bcc/interstitial_00'
+        if not os.path.exists(self.equi_path):
+            os.makedirs(self.equi_path)
+        if not os.path.exists(self.target_path):
+            os.makedirs(self.target_path)
+
+        self.confs = _jdata["structures"]
+        self.inter_param = _jdata["interaction"]
+        self.prop_param = _jdata['properties']
+
+        self.interstitial = Interstitial(_jdata['properties'][0])
+
+    def tearDown(self):
+        if os.path.exists(self.equi_path):
+            shutil.rmtree(self.equi_path)
+        if os.path.exists(self.target_path):
+            shutil.rmtree(self.target_path)
+
+    def test_task_type(self):
+        self.assertEqual('interstitial', self.interstitial.task_type())
+
+    def test_task_param(self):
+        self.assertEqual(self.prop_param[0], self.interstitial.task_param())
+
+    def test_make_confs_bcc(self):
+        if not os.path.exists(os.path.join(self.equi_path, 'CONTCAR')):
+            with self.assertRaises(RuntimeError):
+                self.interstitial.make_confs(self.target_path, self.equi_path)
+        shutil.copy(os.path.join(self.source_path, 'CONTCAR_V_bcc'), os.path.join(self.equi_path, 'CONTCAR'))
+        task_list = self.interstitial.make_confs(self.target_path, self.equi_path)
+        dfm_dirs = glob.glob(os.path.join(self.target_path, 'task.*'))
+        self.assertEqual(len(dfm_dirs), 7)
+
+        incar0 = Incar.from_file(os.path.join('vasp_input', 'INCAR.rlx'))
+        incar0['ISIF'] = 3
+
+        self.assertEqual(os.path.realpath(os.path.join(self.equi_path, 'CONTCAR')),
+                         os.path.realpath(os.path.join(self.target_path, 'POSCAR')))
+        ref_st = Structure.from_file(os.path.join(self.target_path, 'POSCAR'))
+        dfm_dirs.sort()
+        for ii in dfm_dirs[:4]:
+            st_file = os.path.join(ii, 'POSCAR')
+            self.assertTrue(os.path.isfile(st_file))
+            st0 = Structure.from_file(st_file)
+            inter_site = st0[-1]
+            inter = pmg_Interstitial(ref_st, inter_site, charge=0.0)
+            st1 = inter.generate_defect_structure(self.prop_param[0]['supercell'])
+            self.assertEqual(st0, st1)
+
+        for ii in dfm_dirs[4:]:
+            st_file = os.path.join(ii, 'POSCAR')
+            self.assertTrue(os.path.isfile(st_file))
+            st0 = Structure.from_file(st_file)
+            inter_site1 = st0.pop(-1)
+            inter_site2 = st0.pop(-1)
+            center = (inter_site1.coords + inter_site2.coords) / 2
+            self.assertTrue((center[0] - center[1]) < 1e-4)
+            self.assertTrue((center[1] - center[2]) < 1e-4)