Python verlet.VelocityVerlet类代码示例

def maketraj(atoms, t, nstep):
    e = [atoms.get_potential_energy()]
    print "Shape of force:", atoms.get_forces().shape
    dyn = VelocityVerlet(atoms, 5*units.fs)
    for i in range(nstep):
        dyn.run(10)
        energy = atoms.get_potential_energy()
        e.append(energy)
        if ismaster:
            print "Energy: ", energy
        if t is not None:
            t.write()
    return e

        def test_apply_strain(self):
            calc = TersoffScr(**Tersoff_PRB_39_5566_Si_C__Scr)
            timestep = 1.0*units.fs

            atoms = ase.io.read('cryst_rot_mod.xyz')
            atoms.set_calculator(calc)

            # constraints
            top = atoms.positions[:, 1].max()
            bottom = atoms.positions[:, 1].min()
            fixed_mask = ((abs(atoms.positions[:, 1] - top) < 1.0) |
                          (abs(atoms.positions[:, 1] - bottom) < 1.0))
            fix_atoms = FixAtoms(mask=fixed_mask)

            # strain
            orig_height = (atoms.positions[:, 1].max() - atoms.positions[:, 1].min())
            delta_strain = timestep*1e-5*(1/units.fs)
            rigid_constraints = False
            strain_atoms = ConstantStrainRate(orig_height, delta_strain)
            atoms.set_constraint(fix_atoms)

            # dynamics
            np.random.seed(0)
            simulation_temperature = 300*units.kB
            MaxwellBoltzmannDistribution(atoms, 2.0*simulation_temperature)
            dynamics = VelocityVerlet(atoms, timestep)

            def apply_strain(atoms, ConstantStrainRate, rigid_constraints):
                ConstantStrainRate.apply_strain(atoms, rigid_constraints)

            dynamics.attach(apply_strain, 1, atoms, strain_atoms, rigid_constraints)
            dynamics.run(100)

            # tests
            if rigid_constraints == True:
                answer = 0
                temp_answer = 238.2066417638124
            else:
                answer = 0.013228150080099255
                temp_answer = 236.76904696481486

            newpos = atoms.get_positions()
            current_height = newpos[:, 1].max() - newpos[:, 1].min()
            diff_height = (current_height - orig_height)
            self.assertAlmostEqual(diff_height, answer)

            temperature = (atoms.get_kinetic_energy()/(1.5*units.kB*len(atoms)))
            self.assertAlmostEqual(temperature, temp_answer)

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

    calc = CP2K(label='test_H2_MD')
    positions = [(0, 0, 0), (0, 0, 0.7245595)]
    atoms = Atoms('HH', positions=positions, calculator=calc)
    atoms.center(vacuum=2.0)

    # Run MD
    MaxwellBoltzmannDistribution(atoms, 0.5 * 300 * units.kB, force_temp=True)
    energy_start = atoms.get_potential_energy() + atoms.get_kinetic_energy()
    dyn = VelocityVerlet(atoms, 0.5 * units.fs)
    #def print_md():
    #    energy = atoms.get_potential_energy() + atoms.get_kinetic_energy()
    #    print("MD total-energy: %.10feV" %  energy)
    #dyn.attach(print_md, interval=1)
    dyn.run(20)

    energy_end = atoms.get_potential_energy() + atoms.get_kinetic_energy()

    assert energy_start - energy_end < 1e-4
    print('passed test "H2_MD"')

    def MD(self):
        """Molecular Dynamic"""
        from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
        from ase import units
        from ase.md import MDLogger
        from ase.io.trajectory import PickleTrajectory
        from ase.md.langevin import Langevin
        from ase.md.verlet import VelocityVerlet

        dyndrivers = {
            'Langevin': Langevin,
            'None': VelocityVerlet,
        }

        useAsap = False

        mol = self.mol
        temperature = self.definedParams['temperature']
        init_temperature = self.definedParams['init_temperature']
        time_step = self.definedParams['time_step']
        nstep = self.definedParams['nstep']
        nprint = self.definedParams['nprint']
        thermostat = self.definedParams['thermostat']
        prop_file = os.path.join(self.definedParams['workdir'],
                                 self.definedParams['output_prefix']+'.out')
        traj_file = os.path.join(self.definedParams['workdir'],
                                 self.definedParams['output_prefix']+'.traj')

        MaxwellBoltzmannDistribution(mol,init_temperature*units.kB)

        if thermostat == 'None':
            dyn = VelocityVerlet(mol, time_step*units.fs)
        elif thermostat == 'Langevin':
            dyn = Langevin(mol, time_step*units.fs, temperature*units.kB, 0.01 )
        else:
            raise ImplementationError(method,'Thermostat is not implemented in the MD function')

        #Function to print the potential, kinetic and total energy
        traj = PickleTrajectory(traj_file,"a",mol)
        dyn.attach(MDLogger(dyn,mol,prop_file),interval=nprint)
        dyn.attach(traj.write, interval = nprint)

        dyn.run(nstep)
        traj.close()

	else:
		i += 1	
c = FixAtoms(indices = array)
atoms.set_constraint(c)

# relax with Quasi Newtonian
qn = QuasiNewton(atoms, trajectory='qn.traj')
qn.run(fmax=0.001)
write('qn.final.xyz', atoms)

# Set the momenta corresponding to T=300K
MaxwellBoltzmannDistribution(atoms, 300*units.kB)
print 'Removing linear momentum and angular momentum'
Stationary(atoms) # zero linear momentum
ZeroRotation(atoms) # zero angular momentum

# We want to run MD using the VelocityVerlet algorithm.
dyn = VelocityVerlet(atoms, 0.1*units.fs, trajectory='moldyn4.traj') # save trajectory.

#Function to print the potential, kinetic and total energy.
def printenergy(a=atoms):    #store a reference to atoms in the definition.
    epot = a.get_potential_energy() / len(a)
    ekin = a.get_kinetic_energy() / len(a)
    print ("Energy per atom: Epot = %.3feV  Ekin = %.3feV (T=%3.0fK)  Etot = %.3feV" %
           (epot, ekin, ekin/(1.5*units.kB), epot+ekin))
dyn.attach(printenergy, interval=10)

# Now run the dynamics
printenergy()
dyn.run(2000)

    Driver_='VelocityVerlet',
    Driver_MDRestartFrequency=5,
    Driver_Velocities_='',
    Driver_Velocities_empty='<<+ "velocities.txt"',
    Driver_Steps=500,
    Driver_KeepStationary='Yes',
    Driver_TimeStep=8.26,
    Driver_Thermostat_='Berendsen',
    Driver_Thermostat_Temperature=0.00339845142,  # 800 deg Celcius
    # Driver_Thermostat_Temperature=0.0, # 0 deg Kelvin
    Driver_Thermostat_CouplingStrength=0.01)

write_dftb_velocities(test, 'velocities.txt')
os.system('rm md.log.* md.out* geo_end*xyz')
test.set_calculator(calculator_NVE)
dyn = VelocityVerlet(test, 0.000 * fs)  # fs time step.
dyn.attach(MDLogger(dyn, test, 'md.log.NVE', header=True, stress=False,
                    peratom=False, mode='w'), interval=1)
dyn.run(1)  # run NVE ensemble using DFTB's own driver
test = read('geo_end.gen')
write('test.afterNVE.xyz', test)

read_dftb_velocities(test, filename='geo_end.xyz')
write_dftb_velocities(test, 'velocities.txt')

os.system('mv md.out md.out.NVE')
os.system('mv geo_end.xyz geo_end_NVE.xyz')

test.set_calculator(calculator_NVT)
os.system('rm md.log.NVT')
dyn.attach(MDLogger(dyn, test, 'md.log.NVT', header=True, stress=False,

    from asap3 import EMT
    size = 10
else:
    size = 3
    
# Set up a crystal
atoms = FaceCenteredCubic(directions=[[1,0,0],[0,1,0],[0,0,1]], symbol="Cu",
                          size=(size,size,size), pbc=True)

# Describe the interatomic interactions with the Effective Medium Theory
atoms.set_calculator(EMT())

# Set the momenta corresponding to T=300K
MaxwellBoltzmannDistribution(atoms, 300*units.kB)

# We want to run MD with constant energy using the VelocityVerlet algorithm.
dyn = VelocityVerlet(atoms, 5*units.fs)  # 5 fs time step.

#Function to print the potential, kinetic and total energy
def printenergy(a):
    epot = a.get_potential_energy() / len(a)
    ekin = a.get_kinetic_energy() / len(a)
    print ("Energy per atom: Epot = %.3feV  Ekin = %.3feV (T=%3.0fK)  Etot = %.3feV" %
           (epot, ekin, ekin/(1.5*units.kB), epot+ekin))

# Now run the dynamics
printenergy(atoms)
for i in range(20):
    dyn.run(10)
    printenergy(atoms)

    
    qmmm_pot = ForceMixingCarvingCalculator(atoms, qm_region_mask,
                                            mm_pot, qm_pot,
                                            buffer_width=qm_outer_radius,
                                            pbc_type=[False, False, True])
    atoms.set_calculator(qmmm_pot)
#Otherwise it will recover temperature from the previous run.
#Use same random seed so that initialisations are deterministic.
    if not args.restart: 
        print 'Thermalizing atoms'
        np.random.seed(42)
        MaxwellBoltzmannDistribution(atoms, 2.0*sim_T)


    dynamics = VelocityVerlet(atoms, timestep)

    def print_context(ats=atoms, dyn=dynamics):
        print 'steps, T', dyn.nsteps, ats.get_kinetic_energy()/(1.5*units.kB*len(ats))
        print 'G', get_energy_release_rate(ats)/(units.J/units.m**2)
        print 'strain', get_strain(ats)
    dynamics.attach(print_context, interval=8)
    print 'Running Crack Simulation'
    dynamics.run(nsteps)
    print 'Crack Simulation Finished'
  elif args.lotf:
    crack_pos = atoms.info['CrackPos']
    r_scale = 1.00894848312
    mm_pot = Potential('IP EAM_ErcolAd do_rescale_r=T r_scale={0}'.format(r_scale), param_filename=eam_pot, cutoff_skin=2.0)
    #test potential

              (abs(atoms.positions[:, 1] - bottom) < 1.0))
fix_atoms = FixAtoms(mask=fixed_mask)

strain_atoms = ConstantStrainRate(orig_height, strain_rate*timestep)
atoms.set_constraint([fix_atoms, strain_atoms])

mm_pot = Potential("IP SW", cutoff_skin=cutoff_skin)
mm_pot.set_default_properties(['stresses'])

atoms.set_calculator(mm_pot)

# ********* Setup and run MD ***********

MaxwellBoltzmannDistribution(atoms, 2.0*sim_T)

dynamics = VelocityVerlet(atoms, timestep)

def printstatus():
    if dynamics.nsteps == 1:
        print """
State      Time/fs    Temp/K     Strain      G/(J/m^2)  CrackPos/A D(CrackPos)/A 
---------------------------------------------------------------------------------"""

    log_format = ('%(label)-4s%(time)12.1f%(temperature)12.6f'+
                  '%(strain)12.5f%(G)12.4f%(crack_pos_x)12.2f    (%(d_crack_pos_x)+5.2f)')

    atoms.info['label'] = 'D'                  # Label for the status line
    atoms.info['time'] = dynamics.get_time()/units.fs
    atoms.info['temperature'] = (atoms.get_kinetic_energy() /
                                 (1.5*units.kB*len(atoms)))
    atoms.info['strain'] = get_strain(atoms)

                dist = 0.001
            energy += self.A / dist**self.alpha
        return energy

    def __repr__(self):
        return 'Repulsion potential'

    def copy(self):
        return CentralRepulsion(self, R=self.R, A=self.A, alpha=self.alpha)


if __name__ == '__main__':
    from ase.cluster.cubic import FaceCenteredCubic
    from ase.calculators.emt import EMT
    from ase.md.verlet import VelocityVerlet
    from ase.units import fs

    atoms = FaceCenteredCubic('Ag', [(1, 0, 0)], [1], 4.09)
    atoms.center(10)

    atoms.set_calculator(EMT())
    c = ConstantForce(10, [0, 1, 0])  # y=dircted force
    atoms.set_constraint(c)

    md = VelocityVerlet(atoms, 1*fs, trajectory='cf_test.traj',
                        logfile='-')
    md.run(100)

    # from ase.visualize import view
    # view(atoms)

atoms.set_constraint([fix_atoms, strain_atoms])

atoms.set_calculator(params.calc)

# ********* Setup and run MD ***********

# Set the initial temperature to 2*simT: it will then equilibrate to
# simT, by the virial theorem
MaxwellBoltzmannDistribution(atoms, 2.0*params.sim_T)
p = atoms.get_momenta()
p[np.where(fixed_mask)] = 0
atoms.set_momenta(p)

# Initialise the dynamical system
dynamics = VelocityVerlet(atoms, params.timestep)
# dynamics = Langevin(atoms, params.timestep, params.sim_T * units.kB, 1e-4, fixcm=True)

# Print some information every time step
def printstatus():
    if dynamics.nsteps == 1:
        print """
State      Time/fs    Temp/K     Strain      G/(J/m^2)  CrackPos/A D(CrackPos)/A
---------------------------------------------------------------------------------"""

    log_format = ('%(label)-4s%(time)12.1f%(temperature)12.6f'+
                  '%(strain)12.4f%(crack_pos_x)12.2f    (%(d_crack_pos_x)+5.2f)')

    atoms.info['label'] = 'D'                  # Label for the status line
    atoms.info['time'] = dynamics.get_time()/units.fs
    atoms.info['temperature'] = get_temperature(atoms)

      r_scale = 1.00894848312
      pot  = Potential('IP EAM_ErcolAd do_rescale_r=T r_scale={0}'.format(r_scale), param_filename=eam_pot)
      defect.set_calculator(pot)
    else:
      print 'No potential chosen', 1/0

    print 'Finding initial dislocation core positions...'
    try:
      defect.params['core']
    except KeyError:
      defect.params['core'] = np.array([98.0, 98.0, 1.49])

    defect  = set_quantum(defect, params.n_core)
    MaxwellBoltzmannDistribution(defect, 2.0*sim_T)
    if dyn_type =='eam':
       dynamics = VelocityVerlet(defect, timestep)
       dynamics.attach(pass_print_context(defect, dynamics))
    elif dyn_type =='LOTF':
       defect.info['core']= np.array([98.0, 98.0, 1.49])
       print 'Initializing LOTFDynamics'
       verbosity_push(PRINT_VERBOSE)
       dynamics = LOTFDynamics(defect, timestep,
                               params.extrapolate_steps,
                               check_force_error=False)
       dynamics.set_qm_update_func(update_qm_region)
       dynamics.attach(pass_print_context(defect, dynamics))
       dynamics.attach(traj_writer, print_interval, defect)
    else:
      print 'No dyn_type chosen', 1/0
    
    trajectory = AtomsWriter('{0}.traj.xyz'.format(input_file))

def test(temp, frict):
    output = file('Langevin.dat', 'w')
    
    # Make a small perturbation of the momenta
    atoms.set_momenta(1e-6 * np.random.random([len(atoms), 3]))
    print 'Initializing ...'
    predyn = VelocityVerlet(atoms, 0.5)
    predyn.run(2500)

    dyn = Langevin(atoms, timestep, temp, frict)
    print ''
    print ('Testing Langevin dynamics with T = %f eV and lambda = %f' %
           (temp, frict))
    ekin = atoms.get_kinetic_energy()/len(atoms)
    print ekin
    output.write('%.8f\n' % ekin)

    print 'Equilibrating ...'

    # Initial guesses for least-squares fit
    a = 0.04
    b = 2*frict
    c = temp
    params = (a,b,c)
    fitdata = [(0, 2.0 / 3.0 * ekin)]

    tstart = time.time()
    for i in xrange(1,nequil+1):
        dyn.run(nminor)
        ekin = atoms.get_kinetic_energy() / len(atoms)
        fitdata.append((i*nminor*timestep, 2.0/3.0 * ekin))
        if usescipy and i % nequilprint == 0:
            (params, chisq) = leastSquaresFit(targetfunc, params, fitdata)
            print '%.6f  T_inf = %.6f (goal: %f), tau = %.2f,  k = %.6f' % \
                  (ekin, params[2], temp, 1.0/params[1], params[0])
        output.write('%.8f\n' % ekin)
    tequil = time.time() - tstart
    print 'This took %s minutes.' % (tequil / 60)
    output.write('&\n')
    assert abs(temp-params[2]) < 0.25*temp, 'Least-squares fit is way off'
    assert nequil*nminor*timestep > 3.0/params[1], 'Equiliberation was too short'
    fitdata = np.array(fitdata)

    print 'Recording statistical data - this takes ten times longer!'
    temperatures = []
    tstart = time.time()
    for i in xrange(1,nsteps+1):
        dyn.run(nminor)
        ekin = atoms.get_kinetic_energy() / len(atoms)
        temperatures.append(2.0/3.0 * ekin)
        if i % nprint == 0:
            tnow = time.time() - tstart
            tleft = (nsteps-i) * tnow / i
            print '%.6f    (time left: %.1f minutes)' % (ekin, tleft/60)
        output.write('%.8f\n' % ekin)
    output.write('&\n')
    output.close()

    temperatures = np.array(temperatures)
    mean = sum(temperatures) / len(temperatures)
    print 'Mean temperature:', mean, 'eV'
    print
    print 'This test is statistical, and may in rare cases fail due to a'
    print 'statistical fluctuation.'
    print
    assert abs(mean - temp) <= reltol*temp, 'Deviation is too large.'
    print 'Mean temperature:', mean, ' in ', temp, ' +/- ', reltol*temp

    return fitdata, params, temperatures

        print 'No cell geometry given, specifiy either disloc or crack', 1/0 
  
  # ********* Setup and run MD ***********
  # Set the initial temperature to 2*simT: it will then equilibriate to
  # simT, by the virial theorem
  if params.test_mode:
      np.random.seed(0) # use same random seed each time to be deterministic 
  if params.rescale_velo:
      MaxwellBoltzmannDistribution(atoms, 2.0*sim_T)
  # Save frames to the trajectory every `traj_interval` time steps
  # but only when interpolating
  trajectory = AtomsWriter(os.path.join(params.rundir, params.traj_file))
  # Initialise the dynamical system
  system_timer('init_dynamics')
  if params.extrapolate_steps == 1:
      dynamics = VelocityVerlet(atoms, params.timestep)
      check_force_error = False
      if not params.classical:
          qmmm_pot.set(calc_weights=True)
      dynamics.state_label = 'D'
  else:
    print 'Initializing LOTF Dynamics'
    dynamics = LOTFDynamics(atoms, params.timestep,
                            params.extrapolate_steps,
                            check_force_error=check_force_error)
  system_timer('init_dynamics')
#Function to update the QM region at the beginning of each extrapolation cycle   
  if not check_force_error:
      if params.extrapolate_steps == 1:
          if not params.classical:
              dynamics.attach(update_qm_region, 1, dynamics.atoms)

        atoms.set_pbc(False)
        atoms.center(vacuum=6.0)
        atoms.set_velocities(np.zeros_like(atoms.get_positions()))
        cell_c = np.sum(atoms.get_cell()**2, axis=1)**0.5
        N_c = 16 * np.round(cell_c / (0.25 * 16))
        calc = GPAW(gpts=N_c, nbands=1, basis='dzp', setups={'Na': '1'},
                    txt=name + '_gs.txt')
        atoms.set_calculator(calc)
        atoms.get_potential_energy()
        calc.write(name + '_gs.gpw', mode='all')
        del atoms, calc
        time.sleep(10)

    while not os.path.isfile(name + '_gs.gpw'):
        print 'Node %d waiting for file...' % world.rank
        time.sleep(10)
    world.barrier()

    tdcalc = GPAW(name + '_gs.gpw', txt=name + '_td.txt')
    tdcalc.forces.reset() #XXX debug
    tdcalc.initialize_positions()
    atoms = tdcalc.get_atoms()

    traj = PickleTrajectory(name + '_td.traj', 'w', atoms)
    verlet = VelocityVerlet(atoms, timestep * 1e-3 * fs,
                            logfile=paropen(name + '_td.verlet', 'w'),
                            trajectory=traj)
    verlet.attach(Timing(paropen(name + '_td.log', 'w')), ndiv, atoms)
    verlet.run(niter)
    traj.close()