Python VelocityVerlet.run方法代码示例

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import run [as 别名]
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

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import run [as 别名]
        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)

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import run [as 别名]
    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()

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import run [as 别名]
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"')

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import run [as 别名]
	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)

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import run [as 别名]
    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)

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import run [as 别名]
#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

    #quippy using atomic units
    atoms.cutoff = 3.0
    atoms.set_cutoff(3.0)          
    atoms.calc_connect()
    
    qmmm_pot = ForceMixingPotential(pot1=mm_pot, pot2=qm_pot, atoms=atoms,
                                    qm_args_str='carve_cluster=T cluster_periodic_z=T cluster_calc_connect=T '+
                                                'single_cluster=T cluster_vacuum=5.0 terminate=F print_clusters=F '+

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import run [as 别名]
                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)

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import run [as 别名]
    print log_format % atoms.info


dynamics.attach(printstatus)

# Check if the crack has advanced, and stop incrementing the strain if it has
def check_if_cracked(atoms):
    #crack_pos = find_tip_stress_field(atoms)
    # FIXME TS has no local virial
    crack_pos = [0.0, 0.0, 0.0]

    # stop straining if crack has advanced more than tip_move_tol
    if (not atoms.info['is_cracked'] and
        (crack_pos[0] - orig_crack_pos[0]) > params.tip_move_tol):
        atoms.info['is_cracked'] = True
        del atoms.constraints[atoms.constraints.index(strain_atoms)]


dynamics.attach(check_if_cracked, 1, atoms)

# Save frames to the trajectory every `traj_interval` time steps
traj_file = open('traj_verlet_6R.xyz', 'w')
trajectory_write = lambda: ase.io.extxyz.write_xyz(traj_file, atoms)
dynamics.attach(trajectory_write, params.traj_interval)

try:
    # Start running!
    dynamics.run(params.nsteps)
finally:
    traj_file.close()

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import run [as 别名]
    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,
                    peratom=False, mode='w'), interval=1)
dyn.run(1)  # run NVT ensemble using DFTB's own driver
test = read('geo_end.gen')

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import run [as 别名]
        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()

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import run [as 别名]
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(t=atoms):    #store a reference to atoms in the definition.
    epot = t.get_potential_energy() / len(t)
    ekin = t.get_kinetic_energy() / len(t)

    print ("Energy per atom: Epot = %.3feV  Ekin = %.3feV (T=%3.0fK)  Etot = %.3feV" %
           (epot, ekin, ekin/(1.5*units.kB), epot+ekin))
    ETOT = ekin + epot
    print ETOT
    x_dist = atoms.get_distance(a0 = 0, a1 = 1)
    forces = atoms.get_forces()
    OP.append(x_dist)
    OP.append(ekin+epot)
    
dyn.attach(printenergy, interval=10)

# Now run the dynamics
dyn.run(400)
OP_f = [OP[i:i+2] for i in xrange(0,len(OP),2)]
OP_f_text = np.savetxt('out.tmp', OP_f)         # reads array and saves into a file as a string
OP_f_infile = open('out.tmp','r')       # open file to write string array onto
OP_f_text = OP_f_infile.read()

text_file.write(OP_f_text)
text_file.close()

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import run [as 别名]
    def update_atoms(self):
        pass

    def set_atoms(self, atoms):
        self.atoms = atoms
        self.update_atoms()


if __name__ == '__main__':
    from ase import Atoms
    atoms = Atoms(symbols='CeO', cell=[2, 2, 5], positions=np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.4935832]]))
    atoms[0].charge = +4
    atoms[1].charge = -2
    calc = Buck({('Ce', 'O'): (1176.3, 0.381, 0.0)})
    #calc = Buck( {('Ce', 'O'): (0.0, 0.149, 0.0)} )
    atoms.set_calculator(calc)
    print('Epot = ', atoms.get_potential_energy())
    print('Force = ', atoms.get_forces())

    from ase.md.verlet import VelocityVerlet
    dyn = VelocityVerlet(atoms, dt=0.1*units.fs, trajectory='test.traj', logfile='-')
    dyn.run(1000)
    print('coodrs = ', atoms.get_positions())
    from ase.visualize import view
    view(atoms)

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import run [as 别名]
    #    pass

    def __repr__(self):
        return 'Push atoms out of the cell back to the cell'

    def copy(self):
        return ConstantForce(a=self.index, force=self.force)

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
    from constantforce import ConstantForce

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

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

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

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

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import run [as 别名]
    ReportTest("Initial potential energy", epot, -301358.3, 0.5)
    etot = epot + atoms.get_kinetic_energy()

    if 0:
        if cpulayout:
            traj = ParallelNetCDFTrajectory("parallel.nc", atoms)
        else:
            traj = NetCDFTrajectory("serial.nc", atoms)
        traj.Add("PotentialEnergies")
        traj.Update()
        traj.Close()
        print "Trajectory done"
    
    dyn = VelocityVerlet(atoms, 3*units.fs)
    etot2 = None
    for i in range(5):
        dyn.run(15)
        newetot = atoms.get_potential_energy()+ atoms.get_kinetic_energy()
        print >>stdout, "Total energy:", newetot
        temp = atoms.get_kinetic_energy() / (1.5*units.kB*natoms)
        print >>stdout, "Temp:", temp, "K"
        if etot2 == None:
            ReportTest("Total energy (first step)", newetot, etot, 40.0)
            etot2=newetot
        else:
            ReportTest(("Total energy (step %d)" % (i+1,)),
                       newetot, etot2, 20.0)
    print >>stdout, " *** This test completed ***"

ReportTest.Summary()