CN103279645B - Carbon nanotube molecule dynamic-simulation method based on GPU concurrent operation - Google Patents

Abstract

The invention discloses a kind of carbon nanotube molecule dynamic-simulation method based on GPU concurrent operation, comprise the following steps: 1) constructed by CNT and the molecular CNT system model of C60 by database file; 2) C60 molecule is stressed, and vibration back and forth on equilbrium position, when C60 molecule is during in poised state, changes carbon particle on CNT tube wall and position and the speed of C60 molecule; 3) CNT is divided into multi tiered computing unit on CUDA platform, adopts CPU traversal computing unit to carry out trocheameter and calculate that obtain can concurrent operation queue, the stream processing unit of scheduling GPU carries out concurrent operation and processing; 4) repeated execution of steps 3), finally export the running orbit of C60 molecule along with simulated time, and describe the temperature changing curve diagram transmitting along with energy in CNT, complete simulation process. The present invention greatly improves the operation efficiency of Molecular Dynamics.

Description

Carbon nanotube molecule dynamic-simulation method based on GPU concurrent operation

Technical field

The present invention relates to a kind of CNT analytic dynamics emulation mode, especially a kind of based on GPU concurrent operationCarbon nanotube molecule dynamic-simulation method, belongs to simulation technical field.

Background technology

In carbon nanotube molecule dynamics simulation process, owing to relating to ultra-large molecular weight and a large amount of in its computingComplicated floating-point operation, make carbon nanotube molecule dynamics simulation need extremely strong computing capability. Because CNT byA large amount of molecular compositions, relate to a large amount of floating-point operations, along with the expansion of Molecular Dynamics scale, emulation in its analog simulationResult visual also extremely important. For these carbon nanotube molecule dynamics simulation systems, raising simulation algorithmFloating-point operation ability and computation capability have significant role.

In emerging nanometer engineering field, the Macro Mechanism being based upon on continuous media basis is difficult to explain in nanometer engineeringSome phenomena that occur, thereby Molecular Dynamics method becomes one of important research means. Carbon nanotube molecule powerLearn very, the motion of mainly simulating CNT system by Newtonian mechanics, extracts in the data acquisition system in the time of its different conditionsData sample, thereby the configuration integration of calculating CNT system, and further calculate carbon as basis taking the result of configuration integrationThe thermodynamic quantity of nanotube system and other macroscopic properties, thereby efficient emulation CNT system. But carbon nanotube molecule is movingMechanics Simulation algorithm is limited to greatly by computing capability two: in first little space, comprise huge amount particle, amount of calculation is very big; They are two years oldThat typical molecular dynamics time step is femtosecond (fs) level, to computing power requirement in order to ensure numerical simulation stabilityHigh.

Be limited to the development of computing power, for expanding simulation algorithm scale, domestic and international many experts and scholars are to Molecule MotionMechanics Simulation algorithm has done large quantity research. By the serial algorithm that expands at first unit scale, develop into now and appoint by calculatingBusiness is distributed to multiple CPU and is expanded the parallel algorithm of simulation scale, and the scale that can simulate is improved by several thousand initial atomsTo up to a million, the up to ten million scales that even arrive more than one hundred million atoms. Although the computing capability of CPU is very powerful, for toolThe ultra-large simulation system being of practical significance, its speed of service also far can not practical requirement, and for itBatch process carbon, the serial implementation efficiency of single CPU is also lower comparatively speaking. Therefore, large-scale carbon nanotube molecule powerLearning simulation algorithm needs further to break through.

Summary of the invention

The object of the invention is the defect in order to solve above-mentioned prior art, a kind of carbon based on GPU concurrent operation is providedNanotube Molecular Dynamics method, the method can increase substantially the simulation scale of CNT, greatly improves moleculeThe operation efficiency of dynamics simulation.

Object of the present invention can be by taking following technical scheme to reach:

Carbon nanotube molecule dynamic-simulation method based on GPU concurrent operation, is characterized in that comprising the following steps:

1), by database file structure CNT system model, this model is by CNT and be arranged in CNTThe C60 molecular composition of portion, initializes model, the speed of all particles and position coordinates in reading model;

2) C60 molecule is subject to the power that on CNT tube wall, carbon particle produces, and on equilbrium position, vibration back and forth, works as C60Molecule is in the time of poised state, and the condition that different temperatures is set respectively at the forward and backward two ends of CNT is simultaneously simulated, andSet simulated time, speed and the position of carbon particle and C60 molecule on CNT tube wall constantly changed, and then makeStressed, motion state and the energy of C60 molecule constantly change;

3) on CUDA platform, CNT is cut apart, be divided into multilayer computing unit independently mutually, adopt CPUTraversal computing unit carries out trocheameter and calculates that obtain can concurrent operation queue, then dispatches the stream processing unit of GPU, adopts Verlet to calculateMethod is carried out concurrent operation and processing, repeats this step until simulated time finishes;

4) adopt database file to record related data, export the running orbit of C60 molecule along with simulated time, and describeThe temperature changing curve diagram transmitting along with energy in CNT, completes simulation process.

As a kind of preferred version, described CNT is stratiform hollow structure, is made up of multiple hexagonal carbon ring structures;Described C60 molecule has one or more.

As a kind of preferred version, described step 2) in, employing Gaussian heating bath simulation is forward and backward CNTTwo ends arrange respectively the heating bath condition of different temperatures, and set the heating bath time as simulated time.

As a kind of preferred version, described step 3) in, adopt CPU traversal computing unit carry out trocheameter calculator body asUnder:

If a) all particles, all as center calculation particle, jump to step h);

B) find in order first non-competing particle also not by the particle of center calculation, adding can concurrent operation teamRow;

C) all neighbours of this particle of mark are that this parallel queue once can not parallel particle, if neighbour be two degree notCan parallel particle, be once promoted to and can not calculate;

D) all time neighbour of this particle of mark are that this parallel queue two spends can not parallel particle, if inferior neighbour was for once walking abreastParticle, does not revise its number of degrees;

E) this particle of mark is for pressing center calculation particle;

F) judge whether to traverse particle tail of the queue, if continue to carry out, jump to step a) if not;

G) next one can parallel queue start, and returns to step b);

H) finish.

As a kind of preferred version, described step 3) in, adopt Verlet algorithm to carry out concurrent operation and process concrete asUnder:

A) press the position of all particles of CNT system model, calculate the key relation between neighbour's particle and inferior neighbourAnd angular relationship;

B) stream processing unit of scheduling GPU, parallel computation is segmented in the particle in different computing units, and accumulation is calculatedThe interaction force of each particle and its neighbour's particle on CNT tube wall;

C) according to C60 molecule region, calculate in C60 molecule in each particle and its CNT region, place oftenThe interaction force of individual particle;

D), according to the suffered power of particle and speed thereof, upgrade particle position;

E) according to the model moral of particle Interaction Force and the each particle of C60 and tube wall in new positional value calculating tube wallHua Li;

F), according to the suffered power of particle, calculate the heat flow value of this speed of particle and CNT forward and backward two ends heating bath;

If g) reach circulation frequency, calculate and finish; Otherwise the data of preserving particle at CPU end interval, return to stepd)。

As a kind of preferred version, described step 3) in, Verlet algorithm specific design is as follows:

A) by x (t+ Δ t) and x (t-Δ t) carries out Taylor expansion and is shown below:

$\stackrel{\→}{x}(t+\mathrm{\Δ}t)=\stackrel{\→}{x}\left(t\right)+\stackrel{\→}{v}\left(t\right)\mathrm{\Δ}t+\frac{\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2}{2}+\frac{\stackrel{\→}{b}\left(t\right){\mathrm{\Δt}}^3}{2}+O\left({\mathrm{\Δt}}^4\right)---\left(1\right)$

$\stackrel{\→}{x}(t-\mathrm{\Δ}t)=\stackrel{\→}{x}\left(t\right)-\stackrel{\→}{v}\left(t\right)\mathrm{\Δ}t+\frac{\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2}{2}-\frac{\stackrel{\→}{b}\left(t\right){\mathrm{\Δt}}^3}{2}+O\left({\mathrm{\Δt}}^4\right)---\left(2\right)$

Wherein, x (t) is the position of current time t, and Δ t is time step, x (t+ Δ t) represent after the position in a moment, x(t-Δ t) is expressed as the position of previous moment;

B) formula (1) and formula (2) are added, and it is as follows after arranging, to obtain L-expression:

$\stackrel{\→}{x}(t+\mathrm{\Δ}t)=2\stackrel{\→}{x}\left(t\right)-\stackrel{\→}{x}(t-\mathrm{\Δ}t)+\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2+O\left({\mathrm{\Δt}}^4\right)---\left(3\right)$

At position and the acceleration of known CNT system model particle current time t, and previous moment t-Δ tIn the situation of position, the position of a moment t+ Δ t after extrapolating;

C) formula (1) and formula (2) are subtracted each other, then both sides are simultaneously divided by 2 Δ t, the expression formula that obtains speed is as follows:

$\stackrel{\→}{v}\left(t\right)=\frac{\stackrel{\→}{x}(t+\mathrm{\Δ}t)-\stackrel{\→}{x}(t-\mathrm{\Δ}t)}{2\mathrm{\Δ}t}+O\left({\mathrm{\Δt}}^2\right)---\left(4\right)$

In the position of known CNT system model particle previous moment t-Δ t and the position of a rear moment t+ Δ tIn situation, extrapolate the speed v (t) of current time t;

D) by formula (1)～(4), in position, the position in t-Δ t moment and the t-Δ t moment in known particle t-2 Δ t momentThe situation of acceleration under, start Verlet algorithm and carry out integration: according to the t-2 Δ position in t moment, the position in t-Δ t momentWith the t-Δ acceleration in t moment, by t=t-Δ t substitution formula (3), obtain the position of current time t; According to the position of current time tPut, upgrade the acceleration of current time t based on certain potential function; Meanwhile, according to the position of current time t and t-2 Δ t momentPosition, by t=t-Δ t substitution formula (4), upgrade the speed in t-Δ t moment; Obtain position, the t-Δ of particle current time tThe acceleration in the speed in t moment and t moment, repeats this step.

The present invention has following beneficial effect with respect to prior art:

1, the present invention adopts database file structure Nano carbon balls C60 model and carbon nano tube structure model, thereby to carbonNanotube molecular dynamics is carried out emulation, and not only convenient calculating, is also convenient to construct multiple different model and compares, portabilityHigher, and utilize C60 molecule that algorithm can the simulate football shaped structure running orbit in CNT.

2, the present invention adopts the concurrent operation platform CUDA based on GPU being issued by NVIDIA company to carry out image placeReason, CUDA platform is an architecture based on C language, possesses and has a large amount of high-performance calculation instructions and good programming to connectMouthful, can greatly improve the efficiency of carbon nanotube molecule dynamic-simulation method.

3, the present invention can be cut apart large-scale CNT system model by CUDA platform, is divided into multilayer largeLittle suitable, separate and computing unit that GPU can bear, makes the parallel running in CUDA of each algorithm, has greatly improved carbonThe speed of nanotube Molecular Dynamics method.

Brief description of the drawings

Fig. 1 is C60 molecular schematic diagram of the present invention.

Fig. 2 is CNT schematic diagram of the present invention.

Fig. 3 is the structural representation of the present invention's CNT system model of constructing.

Fig. 4 is the data flow block diagram of the present invention's carbon nanotube molecule dynamics simulation in CPU and GPU.

Fig. 5 is the schematic flow sheet that the present invention carries out carbon nanotube molecule dynamics simulation under CUDA environment.

Fig. 6 is comparison diagram running time that the present invention is based on GPU concurrent operation and traditional serial arithmetic.

Detailed description of the invention

Embodiment 1:

As shown in Fig. 1～Fig. 5, the carbon nanotube molecule dynamic-simulation method of the present embodiment based on GPU concurrent operation asUnder:

1) in database file, read the Parameter Conditions of CNT system, as initial temperature, population, density timeDeng, structure CNT system model MOD1 and parameter model MOD2, model M OD2 can arrange various parameters at model M OD1Upper use, model M OD1 as shown in FIG. 1 to 3, is divided with a football shaped C60 who is arranged on CNT inside by CNTSon (molecule being made up of 60 carbon atoms, likeness in form football, has another name called football alkene) composition, described CNT is stratiform hollow knotStructure, pipe shaft is director circle tubular construction, is made up of multiple hexagonal carbon ring structures unit (being made up of carbon atom), its diameter is generally oneBetween tens nanometers, length is far longer than its diameter, this model is initialized to all particles in reading modelSpeed and position coordinates;

2) C60 molecule is subject to the power that on CNT tube wall, carbon particle produces, and on equilbrium position, vibration back and forth, works as C60Molecule, in the time of poised state, adopts Gaussian hot bath method (constraint temperature control method), in the same time-division of CNT two endsThe heating bath condition of different temperatures is not set, and sets the heating bath time, make on CNT tube wall carbon particle and C60 molecule because of balanceState is broken and the variation of occurrence positions and speed, and these variations are interactional, has transitivity, is accompanied by thisThe variation of a little particle positions and speed, energy transmits in the CNT system of elongated tubular, because C60 molecule is positioned at carbonIn nanotube, its stressed meeting, along with the position of CNT constituent particle changes and changes, is therefore followed the mistake of this energy transmissionJourney, the motion state of C60 molecule and energy also can change;

Adopting the general principle of Gaussian hot bath method is in the equation of motion, to add frictional force f_i, and by itself and particle speedDegree v_iConnect, its force bearing formulae is:

f_i＝ma_i+ξmv_i

In the time of equilibrium state, system temperatureConstant, therefore have dEk/dt=0 to have:

Can obtain thus:

3) on CUDA platform, CNT is cut apart, be divided into multilayer suitable size, separate calculating listUnit, partition principle is to improve degree of parallelism in the prerequisite of avoiding too much computing repeatedly, the computing unit of cutting apart must be rationally thickRefinement, because computing unit is excessive, parallel not obvious, computing unit is too small, can cause the calculating of too many unnecessary repetition;

As shown in Figure 4, adopt CPU traversal computing unit carry out trocheameter calculate obtain n can concurrent operation queue, as follows:

If a) all particles, all as center calculation particle, jump to step h);

B) find in order first non-competing particle also not by the particle of center calculation, adding can concurrent operation teamRow;

C) all neighbours of this particle of mark are that this parallel queue once can not parallel particle, if neighbour be two degree notCan parallel particle, be once promoted to and can not calculate;

D) all time neighbour of this particle of mark are that this parallel queue two spends can not parallel particle, if inferior neighbour was for once walking abreastParticle, does not revise its number of degrees;

E) this particle of mark is for pressing center calculation particle;

F) judge whether to traverse particle tail of the queue, if continue to carry out, jump to step a) if not;

G) next one can parallel queue start, and returns to step b);

H) finish;

As shown in Figure 4 and Figure 5, the stream processing unit of scheduling GPU, adopts Verlet algorithm to carry out concurrent operation and processing,As follows:

A) press the position of all particles of CNT system model, calculate the key relation between neighbour's particle and inferior neighbourAnd angular relationship;

B) stream processing unit of scheduling GPU, parallel computation is segmented in the particle in different computing units, and accumulation is calculatedThe interaction force of each particle and its neighbour's particle on CNT tube wall;

C) according to C60 molecule region, calculate in C60 molecule in each particle and its CNT region, place oftenThe interaction force (Van der Waals force) of individual particle;

D), according to the suffered power of particle and speed thereof, upgrade particle position;

E) according to the model moral of particle Interaction Force and the each particle of C60 and tube wall in new positional value calculating tube wallHua Li;

F), according to the suffered power of particle, calculate the heat flow value of this speed of particle and CNT forward and backward two ends heating bath;

If g) reach circulation frequency, calculate and finish; Otherwise the data of preserving particle at CPU end interval, return to stepd)。

4) repeated execution of steps 3), until the heating bath time finishes, adopt database file to record related data, output C60Molecule is along with the running orbit of simulated time, and describes the temperature changing curve diagram transmitting along with energy in CNT, completesSimulation process.

Described step 3) in, because interparticle active force effect of distance is larger, must set one apart from critical value, withJudging its position relationship, when interparticle spacing is less than apart from critical value, is neighbour's particle, when interparticle spacing is greater than distanceCritical value, is inferior neighbour, thinks that its mutual molecular force is negligible.

Described step 3) in, Verlet algorithm specific design is as follows:

A) by x (t+ Δ t) and x (t-Δ t) carries out Taylor expansion and is shown below:

$\stackrel{\→}{x}(t+\mathrm{\Δ}t)=\stackrel{\→}{x}\left(t\right)+\stackrel{\→}{v}\left(t\right)\mathrm{\Δ}t+\frac{\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2}{2}+\frac{\stackrel{\→}{b}\left(t\right){\mathrm{\Δt}}^3}{2}+O\left({\mathrm{\Δt}}^4\right)---\left(1\right)$

$\stackrel{\→}{x}(t-\mathrm{\Δ}t)=\stackrel{\→}{x}\left(t\right)-\stackrel{\→}{v}\left(t\right)\mathrm{\Δ}t+\frac{\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2}{2}-\frac{\stackrel{\→}{b}\left(t\right){\mathrm{\Δt}}^3}{2}+O\left({\mathrm{\Δt}}^4\right)---\left(2\right)$

Wherein, x (t) is the position of current time t, and Δ t is time step, x (t+ Δ t) represent after the position in a moment, x(t-Δ t) is expressed as the position of previous moment;

B) formula (1) and formula (2) are added, and it is as follows after arranging, to obtain L-expression:

$\stackrel{\→}{x}(t+\mathrm{\Δ}t)=2\stackrel{\→}{x}\left(t\right)-\stackrel{\→}{x}(t-\mathrm{\Δ}t)+\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2+O\left({\mathrm{\Δt}}^4\right)---\left(3\right)$

At position and the acceleration of known CNT system model particle current time t, and previous moment t-Δ tIn the situation of position, the position of a moment t+ Δ t after extrapolating;

C) formula (1) and formula (2) are subtracted each other, then both sides are simultaneously divided by 2 Δ t, the expression formula that obtains speed is as follows:

$\stackrel{\→}{v}\left(t\right)=\frac{\stackrel{\→}{x}(t+\mathrm{\Δ}t)-\stackrel{\→}{x}(t-\mathrm{\Δ}t)}{2\mathrm{\Δ}t}+O\left({\mathrm{\Δt}}^2\right)---\left(4\right)$

In the position of known CNT system model particle previous moment t-Δ t and the position of a rear moment t+ Δ tIn situation, extrapolate the speed v (t) of current time t;

D) by formula (1)～(4), in position, the position in t-Δ t moment and the t-Δ t moment in known particle t-2 Δ t momentThe situation of acceleration under, start Verlet algorithm and carry out integration: according to the t-2 Δ position in t moment, the position in t-Δ t momentWith the t-Δ acceleration in t moment, by t=t-Δ t substitution formula (3), obtain the position of current time t; According to the position of current time tPut, upgrade the acceleration of current time t based on certain potential function, thereby obtain that it is stressed; Meanwhile, according to current time t'sThe position in position and t-2 Δ t moment, by t=t-Δ t substitution formula (4), upgrades the speed in t-Δ t moment; Obtain particle currentPosition, the speed in t-Δ t moment and the acceleration in t moment (stressed) of moment t, repeat this step.

As shown in Figure 6, can see the GPU concurrent operation that the present invention adopts, its time of implementation will be far smaller than traditionalSerial arithmetic, therefore, its simulation efficiency improves a lot than traditional serial arithmetic.

Embodiment 2:

The main feature of the present embodiment is: described step 1) in, the C60 molecule that is arranged on CNT inside can be manyIndividual, all the other are with embodiment 1.

The above, be only the preferred embodiment of the invention, but protection scope of the present invention is not limited to this, any ripeKnow those skilled in the art in scope disclosed in this invention, according to technical scheme of the present invention and inventive concept thereofBe equal to and replace or change, all belonged to protection scope of the present invention.

Claims (5)

1. the carbon nanotube molecule dynamic-simulation method based on GPU concurrent operation, is characterized in that comprising the following steps:

1) by database file structure CNT system model, this model is by CNT and be arranged on CNT insideC60 molecular composition, initializes model, the speed of all particles and position coordinates in reading model;

2) C60 molecule is subject to the power that on CNT tube wall, carbon particle produces, and vibration back and forth on equilbrium position, when C60 moleculeIn the time of poised state, the condition that different temperatures is set respectively at the forward and backward two ends of CNT is simultaneously simulated, and setsSimulated time, constantly changes speed and the position of carbon particle and C60 molecule on CNT tube wall, and then C60 is dividedStressed, motion state and the energy of son constantly change;

3) on CUDA platform, CNT is cut apart, be divided into multilayer computing unit independently mutually, adopt CPU traversalComputing unit carries out trocheameter and calculates that obtain can concurrent operation queue, then dispatches the stream processing unit of GPU, adopts Verlet algorithm to enterRow concurrent operation and processing;

4) repeated execution of steps 3), until simulated time finishes, adopt database file to record related data, output C60 moleculeAlong with the running orbit of simulated time, and describe the temperature changing curve diagram transmitting along with energy in CNT, complete emulationProcess;

Described step 3) in, it is specific as follows that employing CPU traversal computing unit carries out trocheameter:

If a) all particles, all as center calculation particle, jump to step h);

B) find in order first non-competing particle also not by the particle of center calculation, adding can concurrent operation queue;

C) all neighbours of this particle of mark are that this parallel queue once can not parallel particle, if neighbour has been that two degree can not be alsoRow particle, was once promoted to and can not calculates;

D) all time neighbour of this particle of mark are that this parallel queue two spends can not parallel particle, if inferior neighbour is the grain that once can not walk abreastSon, does not revise its number of degrees;

E) this particle of mark is for pressing center calculation particle;

F) judge whether to traverse particle tail of the queue, if continue to carry out, jump to step a) if not;

G) next one can parallel queue start, and returns to step b);

H) finish.

2. the carbon nanotube molecule dynamic-simulation method based on GPU concurrent operation according to claim 1, its feature existsIn: described CNT is stratiform hollow structure, is made up of multiple hexagonal carbon ring structures; Described C60 molecule has one or manyIndividual.

3. the carbon nanotube molecule dynamic-simulation method based on GPU concurrent operation according to claim 1 and 2, its spyLevy and be: described step 2) in, adopt Gaussian heating bath simulation, at the forward and backward two ends of CNT, not equality of temperature is set respectivelyThe heating bath condition of degree, and set the heating bath time as simulated time.

4. the carbon nanotube molecule dynamic-simulation method based on GPU concurrent operation according to claim 1, its feature existsIn: described step 3) in, adopt Verlet algorithm carry out concurrent operation and process specific as follows:

A) press the position of all particles of CNT system model, calculate key relation and angle between neighbour's particle and inferior neighbourDegree relation;

B) stream processing unit of scheduling GPU, parallel computation is segmented in the particle in different computing units, and accumulation calculating carbon is receivedThe interaction force of each particle and its neighbour's particle on mitron tube wall;

C), according to C60 molecule region, calculate the each grain in the interior each particle of C60 molecule and its CNT region, placeThe interaction force of son;

D), according to the suffered power of particle and speed thereof, upgrade particle position;

E) according to the Van der Waals of particle Interaction Force and the each particle of C60 and tube wall in new positional value calculating tube wallPower;

F), according to the suffered power of particle, calculate the heat flow value of this speed of particle and CNT forward and backward two ends heating bath;

If g) reach circulation frequency, calculate and finish; Otherwise the data of preserving particle at CPU end interval, return to steps d).

5. the carbon nanotube molecule dynamic-simulation method based on GPU concurrent operation according to claim 4, its feature existsIn: described step 3) in, Verlet algorithm specific design is as follows:

A) by x (t+ Δ t) and x (t-Δ t) carries out Taylor expansion and is shown below:

$\stackrel{\→}{x}(t+\mathrm{\Δ}t)=\stackrel{\→}{x}\left(t\right)+\stackrel{\→}{v}\left(t\right)\mathrm{\Δ}t+\frac{\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2}{2}+\frac{\stackrel{\→}{b}\left(t\right){\mathrm{\Δt}}^3}{2}+O\left({\mathrm{\Δt}}^4\right)---\left(1\right)$

$\stackrel{\→}{x}(t-\mathrm{\Δ}t)=\stackrel{\→}{x}\left(t\right)-\stackrel{\→}{v}\left(t\right)\mathrm{\Δ}t+\frac{\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2}{2}-\frac{\stackrel{\→}{b}\left(t\right){\mathrm{\Δt}}^3}{2}+O\left({\mathrm{\Δt}}^4\right)---\left(2\right)$

Wherein, x (t) is the position of current time t, and Δ t is time step, x (t+ Δ t) represent after the position in a moment, x (t-Δ t) is expressed as the position of previous moment;

B) formula (1) and formula (2) are added, and it is as follows after arranging, to obtain L-expression:

$\stackrel{\→}{x}(t+\mathrm{\Δ}t)=2\stackrel{\→}{x}\left(t\right)-\stackrel{\→}{x}(t-\mathrm{\Δ}t)+\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2+O\left({\mathrm{\Δt}}^4\right)---\left(3\right)$

At position and the acceleration of known CNT system model particle current time t, and the position of previous moment t-Δ tSituation under, the position of a moment t+ Δ t after extrapolating;

C) formula (1) and formula (2) are subtracted each other, then both sides are simultaneously divided by 2 Δ t, the expression formula that obtains speed is as follows:

$\stackrel{\→}{v}\left(t\right)=\frac{\stackrel{\→}{x}(t+\mathrm{\Δ}t)-\stackrel{\→}{x}(t-\mathrm{\Δ}t)}{2\mathrm{\Δ}t}+O\left({\mathrm{\Δt}}^2\right)---\left(4\right)$

In the situation of the position of known CNT system model particle previous moment t-Δ t and the position of a rear moment t+ Δ tUnder, extrapolate the speed v (t) of current time t;

D) by formula (1)～(4), in adding of position, the position in t-Δ t moment and the t-Δ t moment in known particle t-2 Δ t momentIn the situation of speed, start Verlet algorithm and carry out integration: according to the t-2 Δ position in t moment, position and the t-in t-Δ t momentThe acceleration in Δ t moment, by t=t-Δ t substitution formula (3), obtains the position of current time t; According to the position of current time t,Calculate the acceleration of current time t; Meanwhile, according to the position of current time t and the position in t-2 Δ t moment, by t=t-Δ t generationEnter formula (4), upgrade the speed in t-Δ t moment; Obtain position, the speed in t-Δ t moment and the t moment of particle current time tAcceleration, repeat this step.