Python optimize.BFGS类代码示例

def get_optimal_h(atoms, bottom, top, natoms, dyn = False):
    
    # This find the optimal h - when the top is sliding:
    
    if not dyn:
        from scipy.optimize import fmin
        pos_init    =   atoms.positions.copy()
        
        zmax        =   np.amax(atoms.positions[bottom][:,2])
        
        def get_epot(z):
            
            new_pos     =   pos_init
            for iat in range(len(atoms)):
                if top[iat]:
                    new_pos[iat][2] = z + zmax
            
            atoms.positions =   new_pos
            e   = atoms.get_potential_energy()/natoms
            print z, e
            return e  
        
        hmin    =   fmin(get_epot, 3.34)
        emin    =   get_epot(hmin)
        
        atoms.positions = pos_init
        print 'optimal height= %.2f and e=%.2f' %(hmin, emin) 
        return emin, hmin
    else:
        dyn         =   BFGS(atoms)
        dyn.run(fmax=0.03)
        e           =   atoms.get_potential_energy()/natoms
        layers      =   find_layers(atoms.positions)[0]
        hmin        =   layers[1] - layers[0]
        return e, hmin

def get_optimal_h(atoms, natoms, dyn = False, show = False):
    
    # This find the optimal h - when the top is sliding:
    
    if not dyn:
        from scipy.optimize import fmin
        pos_init    =   atoms.positions.copy()
        
        def get_epot(z):
            
            new_pos         =   pos_init
            new_pos[:,2]    =   z 
            atoms.positions =   new_pos
            
            e   = atoms.get_potential_energy()/natoms
            #print   z, e
            return  e  
        
        hmin    =   fmin(get_epot, 3.4, disp = 0)
        emin    =   get_epot(hmin)
        
        atoms.positions = pos_init
        
        if show:    print 'optimal height= %.2f and e=%.2f' %(hmin, emin) 
        return emin, hmin
    else:
        dyn         =   BFGS(atoms)
        dyn.run(fmax=0.03)
        e           =   atoms.get_potential_energy()/natoms
        hmin        =   np.average(atoms.positions[:,2])
        return e, hmin

def test_geopt():
    calc = CP2K(label='test_geopt')
    h2 = molecule('H2', calculator=calc)
    h2.center(vacuum=2.0)
    dyn = BFGS(h2)
    dyn.run(fmax=0.05)
    dist = h2.get_distance(0, 1)
    diff = abs((dist - 1.36733746519) / dist)
    assert(diff < 1e-10)
    print('passed test "geopt"')

def run(sigma, atoms):
    calc = CP2K(label = 'molecules/co-relax/sigma{0}'.format(sigma),
          xc='PBE')
    atoms.set_calculator(calc)
    gopt = BFGS(atoms, logfile=None)
    gopt.run(fmax=1e-2)
    e = atoms.get_potential_energy()
    pos = atoms.get_positions()
    d = ((pos[0] - pos[1])**2).sum()**0.5
    print('{0:1.2f}  {1:1.4f}  {2:1.4f}'.format(sigma, e, d))

 def idpp_interpolate(self, traj='idpp.traj', log='idpp.log', fmax=0.1,
                      optimizer=BFGS):
     d1 = self.images[0].get_all_distances()
     d2 = self.images[-1].get_all_distances()
     d = (d2 - d1) / (self.nimages - 1)
     old = []
     for i, image in enumerate(self.images):
         old.append(image.calc)
         image.calc = IDPP(d1 + i * d)
     opt = BFGS(self, trajectory=traj, logfile=log)
     opt.run(fmax=0.1)
     for image, calc in zip(self.images, old):
         image.calc = calc

def run_neb_calculation(cpu):
    images = [PickleTrajectory('H.traj')[-1]]
    for i in range(nimages):
        images.append(images[0].copy())
    images[-1].positions[6, 1] = 2 - images[0].positions[6, 1]
    neb = NEB(images, parallel=True, world=cpu)
    neb.interpolate()

    images[cpu.rank + 1].set_calculator(Calculator())

    dyn = BFGS(neb)
    dyn.run(fmax=fmax)

    if cpu.rank == 1:
        results.append(images[2].get_potential_energy())

def test_pwscf_calculator():
    if not have_ase():
        skip("no ASE found, skipping test")
    elif not have_pwx():
        skip("no pw.x found, skipping test")
    else:
        pseudo_dir = pj(testdir, prefix, 'pseudo')
        print common.backtick("mkdir -pv {p}; cp files/qe_pseudos/*.gz {p}/; \
            gunzip {p}/*".format(p=pseudo_dir))
        at = get_atoms_with_calc_pwscf(pseudo_dir)

        print "scf"
        # trigger calculation here
        forces = at.get_forces()
        etot = at.get_potential_energy()
        stress = at.get_stress(voigt=False) # 3x3
        
        st = io.read_pw_scf(at.calc.label + '.out')
        assert np.allclose(forces, st.forces)
        assert np.allclose(etot, st.etot)
        assert np.allclose(st.stress, -stress * constants.eV_by_Ang3_to_GPa)
        
        # files/ase/pw.scf.out.start is a norm-conserving LDA struct,
        # calculated with pz-vbc.UPF, so the PBE vc-relax will make the cell
        # a bit bigger
        print "vc-relax"
        from ase.optimize import BFGS
        from ase.constraints import UnitCellFilter
        opt = BFGS(UnitCellFilter(at))
        cell = parse.arr2d_from_txt("""
            -1.97281509  0.          1.97281509
             0.          1.97281509  1.97281509
            -1.97281509  1.97281509  0.""")        
        assert np.allclose(cell, at.get_cell())
        opt.run(fmax=0.05) # run only 2 steps
        cell = parse.arr2d_from_txt("""
            -2.01837531  0.          2.01837531
             0.          2.01837531  2.01837531
            -2.01837531  2.01837531  0""")
        assert np.allclose(cell, at.get_cell())

        # at least 1 backup files must exist: pw.*.0 is the SCF run, backed up
        # in the first iter of the vc-relax
        assert os.path.exists(at.calc.infile + '.0')

def relaxBend(bend, left_idxs, right_idxs, edge, bond, mdrelax):
    
    constraints =   []
    constraints.append(FixAtoms(indices = left_idxs))        
    constraints.append(FixAtoms(indices = right_idxs))        
    #twist       =   twistConst_Rod(bend, 1, edge, bond ,F = 20)
    #twist.set_angle(np.pi/3 + 2./180*np.pi)
    #constraints.append(twist)
    
    add_pot =   LJ_potential_smooth(bend, bond)
    constraints.append(add_pot)
    
    calc    =   LAMMPS(parameters=get_lammps_params()) 
    bend.set_calculator(calc)
    # END CALCULATOR
    
    # RELAX
    bend.set_constraint(constraints)
    dyn     =   BFGS(bend, trajectory = mdrelax)
    dyn.run(fmax=0.05)

def test_lammps_calculator():
    if not have_ase():
        skip("no ASE found, skipping test")
    elif not have_lmp():
        skip("no lammps found, skipping test")
    else:
        at = get_atoms_with_calc_lammps()
        at.rattle(stdev=0.001, seed=int(time.time()))
        common.makedirs(at.calc.directory)
        print common.backtick("cp -v utils/lammps/AlN.tersoff {p}/".format(
            p=at.calc.directory))

        print "scf"
        forces = at.get_forces()
        etot = at.get_potential_energy()
        stress = at.get_stress(voigt=False) # 3x3
        
        st = io.read_lammps_md_txt(at.calc.label + '.out')[0]
        assert np.allclose(forces, st.forces)
        assert np.allclose(etot, st.etot)
        assert np.allclose(st.stress, -stress * constants.eV_by_Ang3_to_GPa,
                           atol=1e-10)
        
        print "relax"
        from ase.optimize import BFGS
        opt = BFGS(at, maxstep=0.04)
        opt.run(fmax=0.001, steps=10)
        coords_frac = parse.arr2d_from_txt("""
            3.3333341909920072e-01    6.6666683819841532e-01    4.4325467247779138e-03
            6.6666681184103216e-01    3.3333362368205072e-01    5.0443254824788963e-01
            3.3333341909918301e-01    6.6666683819838046e-01    3.8356759709402671e-01
            6.6666681184101539e-01    3.3333362368201563e-01    8.8356759861713752e-01
            """)
        assert np.allclose(coords_frac, at.get_scaled_positions(), atol=1e-2)

        # at least 1 backup files must exist
        assert os.path.exists(at.calc.infile + '.0')
        assert os.path.exists(at.calc.outfile + '.0')
        assert os.path.exists(at.calc.dumpfile + '.0')
        assert os.path.exists(at.calc.structfile + '.0')

def main():
    if "ASE_CP2K_COMMAND" not in os.environ:
        raise NotAvailable('$ASE_CP2K_COMMAND not defined')

    calc = CP2K(label='test_H2_GOPT')
    atoms = molecule('H2', calculator=calc)
    atoms.center(vacuum=2.0)

    # Run Geo-Opt
    gopt = BFGS(atoms, logfile=None)
    gopt.run(fmax=1e-6)

    # check distance
    dist = atoms.get_distance(0, 1)
    dist_ref = 0.7245595
    assert (dist - dist_ref) / dist_ref < 1e-7

    # check energy
    energy_ref = -30.7025616943
    energy = atoms.get_potential_energy()
    assert (energy - energy_ref) / energy_ref < 1e-10
    print('passed test "H2_GEO_OPT"')

def test(size, R, nk):

    
    atoms_flat  =   get_square_uCell(size)
    
    view(atoms_flat)
    
    # CALCULATOR FLAT
    calc_f      =   Hotbit(SCC=False, kpts=(nk,nk,1),  \
                           txt= path + 'test_consistency/optimization_flat.cal')
    atoms_flat.set_calculator(calc_f)
    
    opt_f       =   BFGS(atoms_flat)
    opt_f.run(fmax = 0.05)
    e_flat      =   atoms_flat.get_potential_energy()
    

    
    
    atoms_c     =   atoms_flat.copy()
    
    L           =   atoms_c.get_cell().diagonal()
     
    atoms_c.set_cell(L) 
    angle       =   L[1]/R
    atoms_c.rotate('y', np.pi/2)
    atoms_c.translate((-atoms_c[0].x, 0, 0) )
    
    for a in atoms_c:
        r0      =   a.position
        phi     =   r0[1]/L[1]*angle
        a.position[0]   =   R*np.cos(phi)
        a.position[1]   =   R*np.sin(phi)
    
    
    atoms_c       =   Atoms(atoms = atoms_c, container = 'Wedge')
    atoms_c.set_container(angle = angle, height = L[0], physical = False, pbcz = True)
    
    if R < 100:
        view(atoms_c.extended_copy((8,1,3)))
    
    # CALCULATOR Cyl
    calc_c      =   Hotbit(SCC=False, kpts=(nk,1, nk), physical_k = False, \
                           txt= path + 'test_consistency/optimization_cyl.cal')
    atoms_c.set_calculator(calc_c)
    
    opt_c       =   BFGS(atoms_c)
    opt_c.run(fmax = 0.05)
    
    
    e_cyl       =   atoms_c.get_potential_energy()

    print 'R = %.2f' %R
    print 'energy flat     = %.6f' %e_flat
    print 'energy cylinder = %.6f' %e_cyl
    print 'energy dif (e_cylinder - eflat)/nAtoms  = %.6f \n' %(-(e_flat - e_cyl)/len(atoms_flat))
    
    return e_flat, e_cyl, len(atoms_flat)

def ase_minimization(indiv, Optimizer):
    """Function to use built in ASE minimizers to minimize atomic positions in structure.
    Input:
        indiv = Individual class object to be optimized
        Optimizer = Optimizer class object with needed parameters
    Output:
        indiv = Optimized Individual class object.
    """
    if 'MU' in Optimizer.debug:
        debug = True
    else:
        debug = False
    cwd1=os.getcwd()
    olammpsmin = Optimizer.lammps_min
    if Optimizer.lammps_min:
        Optimizer.lammps_min = None
    calc2 = setup_calculator(Optimizer)
#     if 'mass' in indiv[0].get_calculator():
#         mass2 = ['1 '+ str(Optimizer.atomlist[0][2])]
#         if len(Optimizer.atomlist) > 1:
#             for i in range(len(Optimizer.atomlist)-1):
#                 mass2.append(str(i+2) + ' ' + str(Optimizer.atomlist[i+1][2]))
#         calc2=LAMMPS(parameters={ 'pair_style' : Optimizer.pair_style, 'pair_coeff' : Optimizer.pair_coeff , 'mass' : mass2 },files=[ Optimizer.pot_file ])
#     else:
#         calc2=LAMMPS(parameters={ 'pair_style' : Optimizer.pair_style, 'pair_coeff' : Optimizer.pair_coeff},files=[ Optimizer.pot_file ])
    if Optimizer.structure==Defect:
        nat=indiv[0].get_number_of_atoms
        sol=Atoms()
        sol.extend(indiv[0])
        sol.extend(indiv.bulko)
        sol.set_calculator(calc2)
        sol.set_cell(indiv.bulko.get_cell())
        sol.set_pbc(True)
        dyn=BFGS(sol)
        dyn.run(fmax=0.001, steps=2500)
        positions=sol[0:nat].get_positions()
        indiv[0].set_positions(positions)
    else:
        atomsdup=indiv[0].copy()
        atomsdup.set_calculator(calc2)
        dyn=BFGS(indiv[0])
        dyn.run(fmax=0.001, steps=2500)
        positions=atomsdup.get_positions()
        indiv[0].set_positions(positions)
    os.chdir(cwd1)
    calc2.clean()
    Optimizer.lammps_min = olammpsmin
    Optimizer.output.write('ASE Minimization mutation performed on individual = '+repr(indiv.index)+'\n')
    muttype='ASEM'
    if indiv.energy==0:
        indiv.history_index=indiv.history_index+'m'+muttype
    else:
        indiv.history_index=repr(indiv.index)+'m'+muttype
    return indiv
		
		

def run_ase_min(totalsol, fmax=0.01, mxstep=1000, fitscheme='totalenfit', STR=''):
    try:
        dyn=BFGS(totalsol)
        dyn.run(fmax=fmax, steps=mxstep)
    except OverflowError:
        STR+='--- Error: Infinite Energy Calculated - Implement Random shake---\n'
        totalsol.rattle(stdev=0.3)
        dyn=BFGS(totalsol)
        dyn.run(fmax=fmax, steps=mxstep)
    except numpy.linalg.linalg.LinAlgError:
        STR+='--- Error: Singular Matrix - Implement Random shake ---\n'
        totalsol.rattle(stdev=0.2)
        dyn=BFGS(totalsol)
        dyn.run(fmax=fmax, steps=mxstep)
    # Get Energy of Minimized Structure
    en=totalsol.get_potential_energy()
    if fitscheme == 'enthalpyfit':
        pressure=totalsol.get_isotropic_pressure(totalsol.get_stress())
    else:
        pressure=0
    volume = totalsol.get_volume()
    energy=en
    return totalsol, energy, pressure, volume, STR

        E.append(dimer.get_potential_energy())
        F.append(dimer.get_forces())

    F = np.array(F)

    # plt.plot(D, E)

    F1 = np.polyval(np.polyder(np.polyfit(D, E, 7)), D)
    F2 = F[:, :3, 0].sum(1)
    error = abs(F1 - F2).max()

    dimer.constraints = FixInternals(
        bonds=[(r, (0, 2)), (r, (1, 2)),
               (r, (3, 5)), (r, (4, 5))],
        angles=[(a, (0, 2, 1)), (a, (3, 5, 4))])
    opt = BFGS(dimer,
               trajectory=calc.name + '.traj', logfile=calc.name + 'd.log')
    opt.run(0.01)

    e0 = dimer.get_potential_energy()
    d0 = dimer.get_distance(2, 5)
    R = dimer.positions
    v1 = R[1] - R[5]
    v2 = R[5] - (R[3] + R[4]) / 2
    a0 = np.arccos(np.dot(v1, v2) /
                   (np.dot(v1, v1) * np.dot(v2, v2))**0.5) / np.pi * 180
    fmt = '{0:>20}: {1:.3f} {2:.3f} {3:.3f} {4:.1f}'
    print(fmt.format(calc.name, -min(E), -e0, d0, a0))
    assert abs(e0 + eexp) < 0.002
    assert abs(d0 - dexp) < 0.006
    assert abs(a0 - aexp) < 2

    nsteps = '0',
    nstfout = '1',
    nstlog = '1',
    nstenergy = '1',
    nstlist = '1',
    ns_type = 'grid',
    pbc = 'xyz',
    rlist = '1.15',
    coulombtype = 'PME-Switch',
    rcoulomb = '0.8',
    vdwtype = 'shift',
    rvdw = '0.8',
    rvdw_switch = '0.75',
    DispCorr = 'Ener')
CALC_MM.generate_topology_and_g96file()
CALC_MM.generate_gromacs_run_file()

CALC_QMMM = AseQmmmManyqm(nqm_regions = 3, 
                          qm_calculators = [CALC_QM1, CALC_QM2, CALC_QM3], 
                          mm_calculator = CALC_MM,
                          link_info = 'byQM')
#                         link_info = 'byFILE') 

SYSTEM = read_gromos('gromacs_qm.g96')
SYSTEM.set_calculator(CALC_QMMM)
DYN = BFGS(SYSTEM)
DYN.run(fmax = 0.05)

print('exiting fine')
LOG_FILE.close()