naromero

from ase.io import readfrom gpaw import GPAW, FermiDiracfrom gpaw.mixer import Mixerfrom glob import globfrom gpaw.utilities import h2gptsfrom gpaw.poisson import PoissonSolverfrom gpaw.hs_operators import MatrixOperatorfrom gpaw.mpi import worldsystem = read('geom.traj')def adjust_to_h(system, h, idiv=16):    if isinstance(h, float):        h = (h, h, h)    cell = system.get_cell().diagonal()    gpts = (cell / h).astype(int) / idiv * idiv + idiv    cell2 = gpts * h    system.set_cell(cell2)    system.center()    # system.set_pbc(1) # NAR: test for pbc    return gptsgpts = adjust_to_h(system, 0.22)parallel = dict(sl_default=(16, 8, 32),               buffer_size=1024,               band=4,               domain=(4,4,8))B = parallel.get('band', 1)def getkwargs(**kwargs):    kwargs1 = dict(nbands=((len(system) * 6) // B) * B,                   gpts=gpts,                   poissonsolver=PoissonSolver(relax='J', eps=1e-10, nn='M'),                   convergence=dict(eigenstates=1e-9,                                    density=1e-4),                   xc='RPBE',                   txt='gpaw.txt',                   maxiter=360,                   mixer=Mixer(0.02, 9, 200.0),                   parallel=parallel,                   occupations=FermiDirac(0.01, maxiter=4000))    kwargs1.update(**kwargs)    return kwargs1            calc = GPAW(**getkwargs())if world.rank == 0:    import os    if os.path.exists('TERMINATE'):        os.remove('TERMINATE')def checkconv():    decision = 0    if world.rank == 0:        if os.path.exists('TERMINATE'):            decision = 1    decision = world.sum(decision)    if decision:        calc.scf.converged = Truecalc.attach(checkconv, 1)def endcalc():    calc.scf.converged = Truecalc.attach(endcalc, 4)from time import timeclass MaxTime:    def __init__(self, nminutes):        self.starttime = time()        self.maxtime = nminutes * 60            def check(self):        decision = 0        if world.rank == 0:            now = time()            if now - self.starttime > self.maxtime:                decision = 1        decision = world.sum(decision)        if decision:            calc.scf.converged = Truecalc.attach(MaxTime(35).check, 1)system.set_calculator(calc)e = system.get_potential_energy()if world.rank == 0:    fd = open('energy.txt', 'w')    print >> fd, 'e =', e    fd.close()#calc.write('restart.hdf5')import pickleimport numpy as npdef savedensity(calc, nt_G, name):    #nt_G = calc.density.nt_sG[0] # XX    gd = calc.wfs.gd    rank = gd.rank    if calc.wfs.world.rank == 0:        f1 = open('parsize.%s.dat' % name, 'w')        #f2 = open('globalshape.dat', 'w')        print >> f1, '%d %d %d' % tuple(gd.parsize_c)        #print >> f2, '%d %d %d' % tuple(gd.N_c)        f1.close()        #f2.close()    if calc.wfs.bd.rank == 0:        # XXX nspins        densityfile = open('%s.rank%d.spin%d.dat' % (name, gd.rank, 0), 'w')        np.save(densityfile, nt_G)        densityfile.close()#savedensity(calc, calc.density.nt_sG[0], 'density')#world.barrier()#savedensity(calc, calc.hamiltonian.vt_sG[0], 'potential')#saveprojections(calc)