1##################################################################
2##  (c) Copyright 2015-  by Jaron T. Krogel                     ##
3##################################################################
4
5
6#====================================================================#
7#  pwscf.py                                                          #
8#    Nexus interface to the PWSCF simulation code.                   #
9#                                                                    #
10#  Content summary:                                                  #
11#    Pwscf                                                           #
12#      Simulation class for PWSCF.                                   #
13#                                                                    #
14#    generate_pwscf                                                  #
15#      User-facing function to create Pwscf simulation objects.      #
16#                                                                    #
17#====================================================================#
18
19
20import os
21from numpy import array
22from generic import obj
23from physical_system import PhysicalSystem
24from simulation import Simulation
25from pwscf_input import PwscfInput,generate_pwscf_input
26from pwscf_analyzer import PwscfAnalyzer
27from execute import execute
28from debug import ci
29
30
31unique_vdw_functionals = [
32    'optb86b-vdw',
33    'vdw-df3', # optB88+vdW
34    'vdw-df',
35    'vdw-df2',
36    'rev-vdw-df2',
37    'vdw-df-c09',
38    'vdw-df2-c09',
39    'rvv10',
40    ]
41repeat_vdw_functionals = [
42    'vdw-df4', # 'optB86b-vdW'
43    ]
44unique_functionals = [
45    'revpbe',
46    'pw86pbe',
47    'b86bpbe',
48    'pbe0',
49    'hse',
50    'gaup',
51    'pbesol',
52    'pbeq2d',
53    'optbk88',
54    'optb86b',
55    'pbe',
56    'wc',
57    'b3lyp',
58    'pbc',
59    'bp',
60    'pw91',
61    'hcth',
62    'olyp',
63    'tpss',
64    'oep',
65    'hf',
66    'blyp',
67    'lda',
68    'sogga',
69    'm06l',
70    'ev93',
71    ]+unique_vdw_functionals
72repeat_functionals = [
73    'q2d', # pbeq2d
74    'pz', # lda
75    ]+repeat_vdw_functionals
76
77vdw_functionals     = set(unique_vdw_functionals+repeat_vdw_functionals)
78allowed_functionals = set(unique_functionals+repeat_functionals)
79
80
81
82class Pwscf(Simulation):
83    input_type = PwscfInput
84    analyzer_type = PwscfAnalyzer
85    generic_identifier = 'pwscf'
86    application = 'pw.x'
87    application_properties = set(['serial','mpi'])
88    application_results    = set(['charge_density','orbitals','structure','restart'])
89
90    supports_restarts = True # supports restartable, but not force restart yet
91
92    vdw_table = None
93
94    @staticmethod
95    def settings(vdw_table=None):
96        # van der Waals family of functional require the vdW table generated by
97        #  generate_vdW_kernel_table.x: specify 'vdw_table' in settings
98        Pwscf.vdw_table = vdw_table
99    #end def settings
100
101    @staticmethod
102    def restore_default_settings():
103        Pwscf.vdw_table = None
104    #end def restore_default_settings
105
106
107    #def propagate_identifier(self):
108    #    self.input.control.prefix = self.identifier
109    ##end def propagate_identifier
110
111
112    def __init__(self,**sim_args):
113        group_atoms   = sim_args.pop('group_atoms',False)
114        sync_from_scf = sim_args.pop('sync_from_scf',True)
115        Simulation.__init__(self,**sim_args)
116        self.sync_from_scf = False
117        calc = self.input.control.get('calculation',None)
118        if calc=='nscf':
119            self.sync_from_scf = sync_from_scf
120        #end if
121        if group_atoms and isinstance(self.system,PhysicalSystem):
122            self.warn('requested grouping by atomic species, but pwscf does not group atoms anymore!')
123            #self.system.structure.group_atoms()
124        #end if
125    #end def __init__
126
127
128    def write_prep(self):
129        #make sure the output directory exists
130        outdir = os.path.join(self.locdir,self.input.control.outdir)
131        if not os.path.exists(outdir):
132            os.makedirs(outdir)
133        #end if
134        #copy over vdw_table for vdW-DF functional
135        if self.path_exists('input/system/input_dft'):
136            functional = self.input.system.input_dft.lower()
137            if '+' not in functional and functional not in allowed_functionals:
138                self.warn('functional "{0}" is unknown to pwscf'.format(functional))
139            #end if
140            if functional in vdw_functionals:
141                if self.vdw_table is None:
142                    self.error('attempting to run vdW functional "{0}", but vdw_table is missing\nplease provide path to table file via "vdw_table" parameter in settings'.format(functional))
143                #end if
144                cd_rel = os.path.relpath(self.vdw_table,self.locdir)
145                # copy instead of link to vdw_table to avoid file-lock from multiple pw.x instances
146                cp_cmd = 'cd '+self.locdir+';cp '+cd_rel+' .'
147                os.system(cp_cmd)
148            #end if
149        #end if
150    #end def write_prep
151
152
153    def check_result(self,result_name,sim):
154        input = self.input
155        control = input.control
156        if result_name=='charge_density' or result_name=='restart':
157            calculating_result = True
158        elif result_name=='orbitals':
159            calculating_result = 'calculation' not in control or 'scf' in control.calculation.lower()
160        elif result_name=='structure':
161            calculating_result = 'calculation' in control and 'relax' in control.calculation.lower()
162        else:
163            calculating_result = False
164        #end if
165        return calculating_result
166    #end def check_result
167
168
169    def get_result(self,result_name,sim):
170        result = obj()
171        input = self.input
172        control = input.control
173        prefix = 'pwscf'
174        outdir = './'
175        if 'prefix' in control:
176            prefix = control.prefix
177        #end if
178        if 'outdir' in control:
179            outdir = control.outdir
180        #end if
181        if outdir.startswith('./'):
182            outdir = outdir[2:]
183        #end if
184        if result_name=='charge_density' or result_name=='restart':
185            result.locdir   = self.locdir
186            result.outdir   = os.path.join(self.locdir,outdir)
187            result.location = os.path.join(self.locdir,outdir,prefix+'.save','charge-density.dat')
188            result.spin_location = os.path.join(self.locdir,outdir,prefix+'.save','spin-polarization.dat')
189        elif result_name=='orbitals':
190            result.location = os.path.join(self.locdir,outdir,prefix+'.wfc1')
191        elif result_name=='structure':
192            pa = self.load_analyzer_image()
193            structs = pa.structures
194            struct  = structs[len(structs)-1]
195            pos   = struct.positions
196            atoms = struct.atoms
197            if 'celldm(1)' in self.input.system:
198                scale = self.input.system['celldm(1)']
199            else:
200                scale = 1.0
201            #end if
202            pos   = scale*array(pos)
203
204            structure = self.system.structure.copy()
205            structure.change_units('B')
206            structure.pos = pos
207            structure.set_elem(atoms)
208            if 'axes' in struct:
209                structure.axes = struct.axes.copy()
210            #end if
211            result.structure = structure
212        else:
213            self.error('ability to get result '+result_name+' has not been implemented')
214        #end if
215        return result
216    #end def get_result
217
218
219    def incorporate_result(self,result_name,result,sim):
220        if result_name=='charge_density':
221            c = self.input.control
222            res_path = os.path.abspath(result.locdir)
223            loc_path = os.path.abspath(self.locdir)
224            if res_path==loc_path:
225                None # don't need to do anything if in same directory
226            elif self.sync_from_scf: # rsync output into nscf dir
227                outdir = os.path.join(self.locdir,c.outdir)
228                command = 'rsync -av {0}/* {1}/'.format(result.outdir,outdir)
229                if not os.path.exists(outdir):
230                    os.makedirs(outdir)
231                #end if
232                sync_record = os.path.join(outdir,'nexus_sync_record')
233                if not os.path.exists(sync_record):
234                    execute(command)
235                    f = open(sync_record,'w')
236                    f.write('\n')
237                    f.close()
238                #end if
239            else: # attempt to use symbolic links instead
240                link_loc = os.path.join(self.locdir,c.outdir,c.prefix+'.save')
241                cd_loc = result.location
242                cd_rel = os.path.relpath(cd_loc,link_loc)
243                sp_loc = result.spin_location
244                sp_rel = os.path.relpath(sp_loc,link_loc)
245                cwd = os.getcwd()
246                if not os.path.exists(link_loc):
247                    os.makedirs(link_loc)
248                #end if
249                os.chdir(link_loc)
250                os.system('ln -s '+cd_rel+' charge-density.dat')
251                os.system('ln -s '+sp_rel+' spin-polarization.dat')
252                os.chdir(cwd)
253            #end if
254        elif result_name=='structure':
255            relstruct = result.structure.copy()
256            relstruct.change_units('B')
257            self.system.structure = relstruct
258            self.system.remove_folded()
259
260            input = self.input
261            preserve_kp = 'k_points' in input and 'specifier' in input.k_points and (input.k_points.specifier=='automatic' or input.k_points.specifier=='gamma')
262            if preserve_kp:
263                kp = input.k_points.copy()
264            #end if
265            input.incorporate_system(self.system)
266            if preserve_kp:
267                input.k_points = kp
268            #end if
269        elif result_name=='restart':
270            c = self.input.control
271            if('startingwfc' in self.input.electrons and self.input.electrons.startingwfc != 'file'):
272                self.error('Exiting. User has specified startingwfc=\''+self.input.electrons.startingwfc+'\'.\nThis value will be overwritten when incorporating result \'restart\'.\nPlease fix conflict.')
273            #end if
274            if('startingpot' in self.input.electrons and self.input.electrons.startingpot != 'file'):
275                self.error('Exiting. User has specified startingpot=\''+self.input.electrons.startingpot+'\'.\nThis value will be overwritten when incorporating result \'restart\'.\nPlease fix conflict.')
276            #end if
277            c.restart_mode='restart'
278            res_path = os.path.abspath(result.locdir)
279            loc_path = os.path.abspath(self.locdir)
280            if res_path==loc_path:
281                None # don't need to do anything if in same directory
282            else: # rsync output into new scf dir
283                outdir = os.path.join(self.locdir,c.outdir)
284                command = 'rsync -av {0}/* {1}/'.format(result.outdir,outdir)
285                if not os.path.exists(outdir):
286                    os.makedirs(outdir)
287                #end if
288                sync_record = os.path.join(outdir,'nexus_sync_record')
289                if not os.path.exists(sync_record):
290                    execute(command)
291                    f = open(sync_record,'w')
292                    f.write('\n')
293                    f.close()
294                #end if
295            #end if
296        else:
297            self.error('ability to incorporate result '+result_name+' has not been implemented')
298        #end if
299    #end def incorporate_result
300
301
302    def check_sim_status(self):
303        outfile = os.path.join(self.locdir,self.outfile)
304        fobj = open(outfile,'r')
305        output = fobj.read()
306        fobj.close()
307        not_converged = 'convergence NOT achieved'  in output
308        time_exceeded = 'Maximum CPU time exceeded' in output
309        user_stop     = 'Program stopped by user request' in output
310        run_finished  = 'JOB DONE' in output
311        restartable = not_converged or time_exceeded or user_stop
312        restart = run_finished and self.restartable and restartable
313        if restart:
314            self.save_attempt()
315            self.input.control.restart_mode = 'restart'
316            self.reset_indicators()
317        else:
318            error_in_routine = 'Error in routine' in output
319            failed = not_converged or time_exceeded or user_stop
320            failed |= error_in_routine
321            self.finished = run_finished
322            self.failed   = failed
323        #end if
324    #end def check_sim_status
325
326
327    def get_output_files(self):
328        output_files = []
329        return output_files
330    #end def get_output_files
331
332
333    def app_command(self):
334        return self.app_name+' -input '+self.infile
335    #end def app_command
336#end class Pwscf
337
338
339
340def generate_pwscf(**kwargs):
341    sim_args,inp_args = Pwscf.separate_inputs(kwargs)
342
343    if not 'input' in sim_args:
344        input_type = inp_args.delete_optional('input_type','generic')
345        sim_args.input = generate_pwscf_input(input_type,**inp_args)
346    #end if
347    pwscf = Pwscf(**sim_args)
348
349    return pwscf
350#end def generate_pwscf
351