Coverage for gpaw/test/big/dcdft/pbe_abinit_hgh.py: 0%

1import os 

2import sys 

3import time 

4 

5import numpy as np 

6 

7import ase.db 

8from ase.units import Rydberg 

9from ase.utils import opencew 

10from ase.calculators.calculator import kpts2mp 

11from ase.io import Trajectory 

12from ase.calculators.abinit import Abinit 

13from ase.test.tasks.dcdft import DeltaCodesDFTCollection as Collection 

14 

15collection = Collection() 

16 

17if len(sys.argv) == 1: 

18 names = collection.names 

19else: 

20 names = [sys.argv[1]] 

21 

22c = ase.db.connect('dcdft_abinit_hgh.db') 

23 

24ecut = 250 

25 

26kptdensity = 16.0 # this is converged 

27kptdensity = 6.0 # just for testing 

28width = 0.01 

29ecutsm = 0.0 

30fband = 1.5 

31tolsym = 1.e-12 

32 

33linspace = (0.98, 1.02, 5) # eos numpy's linspace 

34linspacestr = ''.join([str(t) + 'x' for t in linspace])[:-1] 

35 

36code = 'abinit' + '-' + '_c' + str(ecut) + '_e' + linspacestr 

37code = (code + '_k' + str(kptdensity) + '_w' + str(width) + '_s' + 

38 str(ecutsm) + '_t' + str(tolsym)) 

39 

40for name in names: 

41 # save all steps in one traj file in addition to the database 

42 # we should only used the database c.reserve, but here 

43 # traj file is used as another lock ... 

44 fd = opencew(name + '_' + code + '.traj') 

45 if fd is None: 

46 continue 

47 traj = Trajectory(name + '_' + code + '.traj', 'w') 

48 atoms = collection[name] 

49 if name == 'Mn': # fails to find the right magnetic state 

50 atoms.set_initial_magnetic_moments([10., 20., -10., -20.]) 

51 if name == 'Co': # fails to find the right magnetic state 

52 atoms.set_initial_magnetic_moments([10., 10.]) 

53 if name == 'Ni': # fails to find the right magnetic state 

54 atoms.set_initial_magnetic_moments([10., 10., 10., 10.]) 

55 cell = atoms.get_cell() 

56 kpts = tuple(kpts2mp(atoms, kptdensity, even=True)) 

57 kwargs = {} 

58 # loop over EOS linspace 

59 for n, x in enumerate(np.linspace(linspace[0], linspace[1], linspace[2])): 

60 id = c.reserve(name=name, ecut=ecut, 

61 linspacestr=linspacestr, 

62 kptdensity=kptdensity, width=width, ecutsm=ecutsm, 

63 fband=fband, tolsym=tolsym, 

64 x=x) 

65 if id is None: 

66 continue 

67 # perform EOS step 

68 atoms.set_cell(cell * x, scale_atoms=True) 

69 # set calculator 

70 atoms.calc = Abinit( 

71 pps='hgh.k', # uses highest valence hgh.k pps 

72 label=name + '_' + code + '_' + str(n), 

73 xc='PBE', 

74 kpts=kpts, 

75 ecut=ecut * Rydberg, 

76 occopt=3, 

77 tsmear=width, 

78 ecutsm=ecutsm, 

79 toldfe=1.0e-6, 

80 nstep=900, 

81 pawovlp=-1, # bypass overlap check 

82 fband=fband, 

83 # https://forum.abinit.org/viewtopic.php?f=8&t=35 

84 chksymbreak=0, 

85 tolsym=tolsym, 

86 prtwf=0, 

87 prtden=0, 

88 ) 

89 atoms.calc.set(**kwargs) # remaining calc keywords 

90 t = time.time() 

91 atoms.get_potential_energy() 

92 c.write(atoms, 

93 name=name, ecut=ecut, 

94 linspacestr=linspacestr, 

95 kptdensity=kptdensity, width=width, ecutsm=ecutsm, 

96 fband=fband, tolsym=tolsym, 

97 x=x, 

98 time=time.time() - t) 

99 traj.write(atoms) 

100 wfk = name + '_' + code + '_' + str(n) + 'o_WFK' 

101 if os.path.exists(wfk): 

102 os.remove(wfk) 

103 del c[id]