Python VelocityVerlet.attach方法代码示例

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

    #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,

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import attach [as 别名]
    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)
    atoms.info['G'] = get_energy_release_rate(atoms)/(units.J/units.m**2)
    
    crack_pos = find_crack_tip_stress_field(atoms, calc=mm_pot)
    atoms.info['crack_pos_x'] = crack_pos[0]
    atoms.info['d_crack_pos_x'] = crack_pos[0] - orig_crack_pos[0]

    print log_format % atoms.info
dynamics.attach(printstatus)

def check_if_cracked(atoms):
    crack_pos = find_crack_tip_stress_field(atoms, calc=mm_pot)

    # 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]) > tip_move_tol:
        atoms.info['is_cracked'] = True
        del atoms.constraints[atoms.constraints.index(strain_atoms)]
dynamics.attach(check_if_cracked, 1, atoms)

trajectory = AtomsWriter(traj_file)
dynamics.attach(trajectory, traj_interval, atoms)

dynamics.run(nsteps)

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import attach [as 别名]
    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)
    atoms.info['strain'] = get_strain(atoms)
    # FIXME no local virial in TS, need another way to track the crack
    atoms.info['crack_pos_x'] = 0.
    atoms.info['d_crack_pos_x'] = 0.

    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)

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

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

# 需要导入模块: from ase.md.verlet import VelocityVerlet [as 别名]
# 或者: from ase.md.verlet.VelocityVerlet import attach [as 别名]
      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)
      else:
      # choose appropriate update function for defects or crack.
      # or grainboundary.
        print 'Setting Update Function'
        if geom =='disloc':
          dynamics.set_qm_update_func(update_qm_region)
        elif geom =='crack':
          dynamics.set_qm_update_func(update_qm_region_crack)
        else:
          print 'No geometry chosen', 1/0
  
  if check_force_error:
      pred_corr_logfile = open(os.path.join(params.rundir, 'pred-corr-error.txt'), 'w')
      dynamics.attach(log_pred_corr_errors, 1, dynamics, pred_corr_logfile)
     

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