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

1import sys 

2import time 

3 

4import numpy as np 

5 

6import ase.db 

7from ase.units import Rydberg 

8from ase.utils import opencew 

9from ase.calculators.calculator import kpts2mp 

10from ase.io import Trajectory 

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

12from gpaw import GPAW, PW, FermiDirac 

13from gpaw.utilities import h2gpts 

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_gpaw_pw.db') 

23 

24# mode = 'lcao' 

25mode = 'fd' 

26mode = 'pw' 

27 

28e = 0.08 # h -> gpts 

29e = round(100 * Rydberg, 0) 

30 

31kptdensity = 16.0 # this is converged 

32kptdensity = 6.0 # just for testing 

33width = 0.01 

34 

35relativistic = True 

36constant_basis = True 

37 

38if relativistic: 

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

40else: 

41 linspace = (0.92, 1.08, 7) # eos numpy's linspace 

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

43 

44code = ('gpaw' + '-' + mode + str(e) + '_c' + 

45 str(constant_basis) + '_e' + linspacestr) 

46code = code + '_k' + str(kptdensity) + '_w' + str(width) 

47code = code + '_r' + str(relativistic) 

48 

49for name in names: 

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

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

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

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

54 if fd is None: 

55 continue 

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

57 atoms = collection[name] 

58 cell = atoms.get_cell() 

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

60 kwargs = {} 

61 if mode in ['fd', 'lcao']: 

62 if constant_basis: 

63 # gives more smooth EOS in fd mode 

64 kwargs.update({'gpts': h2gpts(e, cell)}) 

65 else: 

66 kwargs.update({'h': e}) 

67 elif mode == 'pw': 

68 if constant_basis: 

69 kwargs.update({'mode': PW(e, cell=cell)}) 

70 kwargs.update({'gpts': h2gpts(0.10, cell)}) 

71 else: 

72 kwargs.update({'mode': PW(e)}) 

73 if mode == 'pw': 

74 if name in ['Li', 'Na']: 

75 # listserv.fysik.dtu.dk/pipermail/gpaw-developers/2012-May/002870.html 

76 if constant_basis: 

77 kwargs.update({'gpts': h2gpts(0.05, cell)}) 

78 else: 

79 kwargs.update({'h': 0.05}) 

80 if mode == 'lcao': 

81 kwargs.update({'mode': 'lcao'}) 

82 kwargs.update({'basis': 'dzp'}) 

83 if mode == 'fd': 

84 kwargs.update({'mode': 'fd'}) 

85 if name in ['He', 'Ne', 'Ar', 'Kr', 'Xe', 'Rn', 'Ca', 'Sr', 'Ba', 'Be']: 

86 # results wrong / scf slow with minimal basis 

87 kwargs.update({'basis': 'dzp'}) 

88 kwargs.update({'nbands': -5}) 

89 # loop over EOS linspace 

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

91 id = c.reserve(name=name, mode=mode, e=e, linspacestr=linspacestr, 

92 kptdensity=kptdensity, width=width, 

93 relativistic=relativistic, 

94 constant_basis=constant_basis, 

95 x=x) 

96 if id is None: 

97 continue 

98 # perform EOS step 

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

100 # set calculator 

101 base_kwargs = dict( 

102 txt=name + '_' + code + '_' + str(n) + '.txt', 

103 xc='PBE', 

104 kpts=kpts, 

105 occupations=FermiDirac(width), 

106 parallel={'band': 1}, 

107 maxiter=777) 

108 atoms.calc = GPAW(**dict(base_kwargs, **kwargs)) 

109 t = time.time() 

110 atoms.get_potential_energy() 

111 c.write(atoms, 

112 name=name, mode=mode, e=e, linspacestr=linspacestr, 

113 kptdensity=kptdensity, width=width, 

114 relativistic=relativistic, 

115 constant_basis=constant_basis, 

116 x=x, 

117 niter=atoms.calc.get_number_of_iterations(), 

118 time=time.time() - t) 

119 traj.write(atoms) 

120 del c[id]