nano_mpi/src/dynamics_velocity_verlet.F90

module dynamics_velocity_verlet
  use simulation
  use definitions, ONLY : DP, NDIM
  use constants
  use parameter
  use force_module, ONLY : calc_forces
  use velocity,     ONLY : scale_velocities
  use force_utilities, ONLY : check
  use dynamics_utilities, ONLY : evolve_time, init_dynamics_loop, sim_status, &
       DT2,DTSQRT,DTSQ2, init_status
#ifdef MPI
  use mpi_module
#endif

  implicit none


  public  :: velocity_verlet_dynamics
  private :: movea, moveb
 
  type (sim_status) :: verlet_status
  

contains


subroutine velocity_verlet_dynamics()
  use velocity, ONLY: calc_temp

 
  logical                          :: update_nlist
  logical                          :: do_pot
  logical                          :: not_done
  logical                          :: nmflag

  integer                          :: step, i
 
  call init_status(verlet_status)
  call init_dynamics_loop(verlet_status)

  !     Start the main simulation loop.

  not_done = .true.
  nmflag  =  .false.
  step = 0

  
  dynamics: do  
     if (.not.not_done) exit dynamics

     step = step + 1


     if (checktemp.and.(mod(step,check_temp_steps).eq.0)) then
          call calc_temp(verlet_status%temp)
        if (dabs(verlet_status%temp-target_temp).gt.therm_variance) then
           call scale_velocities(target_temp)
        end if
     end if

     call movea()


     call check(update_nlist)


     if (do_pot) then
        call calc_forces(update_nlist,nmflag,pe = verlet_status%pot_e)
     else
        call calc_forces(update_nlist,nmflag)
     endif

     call moveb()         
     

     call evolve_time(step,verlet_status,do_pot,not_done)
 

#ifdef MPI
     call mpi_barrier(mpi_comm_world,mpi_err)
#endif
     end do dynamics
end subroutine velocity_verlet_dynamics

!    *******************************************************************
!    ** FICHE F.4                                                     **
!    ** TWO ROUTINES THAT TOGETHER IMPLEMENT VELOCITY VERLET METHOD.  **
!    **                                                               **
!    ** REFERENCE:                                                    **
!    **                                                               **
!    ** SWOPE ET AL., J. CHEM. PHYS. 76, 637, 1982.                   **
!    **                                                               **
!    ** ROUTINES SUPPLIED:                                            **
!    **                                                               **
!    ** SUBROUTINE MOVEA ( DT, M )                                    **
!    **    MOVES POSITIONS AND PARTIALLY UPDATES VELOCITIES.          **
!    ** SUBROUTINE MOVEB ( DT, M, K )                                 **
!    **    COMPLETES VELOCITY MOVE AND CALCULATES KINETIC ENERGY.     **
!    **                                                               **
!    ** PRINCIPAL VARIABLES:                                          **
!    **                                                               **
!    ** INTEGER N                   NUMBER OF MOLECULES               **
!    ** REAL    DT                  TIMESTEP                          **
!    ** REAL    M                   ATOMIC MASS                       **
!    ** REAL    K                   KINETIC ENERGY                    **
!    ** REAL    RX(N),RY(N),RZ(N)   POSITIONS                         **
!    ** REAL    VX(N),VY(N),VZ(N)   VELOCITIES                        **
!    ** REAL    FX(N),FY(N),FZ(N)   FORCES                            **
!    **                                                               **
!    ** USAGE:                                                        **
!    **                                                               **
!    ** AT THE START OF A TIMESTEP, MOVEA IS CALLED TO ADVANCE THE    **
!    ** POSITIONS AND 'HALF-ADVANCE' THE VELOCITIES.  THEN THE FORCE  **
!    ** ROUTINE IS CALLED, AND THIS IS FOLLOWED BY MOVEB WHICH        **
!    ** COMPLETES THE ADVANCEMENT OF VELOCITIES.                      **
!    *******************************************************************


subroutine movea( )

  !     *******************************************************************
  !     ** FIRST PART OF VELOCITY VERLET ALGORITHM                       **
  !     **                                                               **
  !     ** USAGE:                                                        **
  !     **                                                               **
  !     ** THE FIRST PART OF THE ALGORITHM IS A TAYLOR SERIES WHICH      **
  !     ** ADVANCES POSITIONS FROM T TO T + DT AND VELOCITIES FROM       **
  !     ** T TO T + DT/2.  AFTER THIS, THE FORCE ROUTINE IS CALLED.      **
  !     *******************************************************************


  integer ::     i, dim
  !
  !     units for time are femtosec  (10^-15 sec)
  !     units for distance are angstroms (10^-10 m)
  !     units for velocity are angstroms femtosec^-1
  !     units for mass are a.m.u. (1.661 * 10^-27 kg)
  !     units for force are kcal mol^-1 angstrom^-1
  !
  !     converter will put the final terms into angstroms.
  !       or angstrom/femtosecond.

!  converter = 1.0d0/2.390664d3

  !     *******************************************************************


  do i = 1,nlocal
     do dim = 1, 3
        
        q(dim,i) = q(dim,i) + dt*v(dim,i) + (dtsq2*f(dim,i)/mass(i))*KCALMOL_TO_AMUA2FS2
        
        v(dim,i) = v(dim,i) + (dt2*f(dim,i)/mass(i))*KCALMOL_TO_AMUA2FS2
        
     end do
  end do
  
end subroutine movea
      
      
subroutine moveb()
#ifdef MPI
  use mpi_module
#endif
  !     *******************************************************************
  !     ** SECOND PART OF VELOCITY VERLET ALGORITHM                      **
  !     **                                                               **
  !     ** USAGE:                                                        **
  !     **                                                               **
  !     ** THE SECOND PART OF THE ALGORITHM ADVANCES VELOCITIES FROM     **
  !     ** T + DT/2 TO T + DT. THIS ASSUMES THAT FORCES HAVE BEEN        **
  !     ** COMPUTED IN THE FORCE ROUTINE AND STORED IN FX, FY, FZ.       **
  !     *******************************************************************


  integer :: i, dim
  !
  !     units for time are femtosec  (10^-15 sec)
  !     units for distance are angstroms (10^-10 m)
  !     units for velocity are angstroms femtosec^-1
  !     units for mass are a.m.u. (1.661 * 10^-27 kg)
  !     units for force are kcal mol^-1 angstrom^-1
  !
  !     converter will put the final terms into kcal/mol
  !       or angstrom/femtosecond.

!  converter = 1.0d0/2.390664d3

  !     *******************************************************************

 
  do i = 1, nlocal
     do dim = 1, 3
        v(dim,i) = v(dim,i) + (dt2*f(dim,i)/mass(i))*KCALMOL_TO_AMUA2FS2
     end do
  end do

end subroutine moveb

!! modified to be used in a mpi simulation.....


end module dynamics_velocity_verlet
Revision:	4
Committed:	Mon Jun 10 17:18:36 2002 UTC (22 years, 1 month ago) by chuckv
File size:	7102 byte(s)
Log Message:	Import Root
#	User	Rev	Content
1	chuckv	4	module dynamics_velocity_verlet
2			use simulation
3			use definitions, ONLY : DP, NDIM
4			use constants
5			use parameter
6			use force_module, ONLY : calc_forces
7			use velocity, ONLY : scale_velocities
8			use force_utilities, ONLY : check
9			use dynamics_utilities, ONLY : evolve_time, init_dynamics_loop, sim_status, &
10			DT2,DTSQRT,DTSQ2, init_status
11			#ifdef MPI
12			use mpi_module
13			#endif
14
15			implicit none
16
17
18
19
20			public :: velocity_verlet_dynamics
21			private :: movea, moveb
22
23			type (sim_status) :: verlet_status
24
25
26
27			contains
28
29
30
31
32			subroutine velocity_verlet_dynamics()
33			use velocity, ONLY: calc_temp
34
35
36
37			logical :: update_nlist
38			logical :: do_pot
39			logical :: not_done
40			logical :: nmflag
41
42			integer :: step, i
43
44			call init_status(verlet_status)
45			call init_dynamics_loop(verlet_status)
46
47			! Start the main simulation loop.
48
49			not_done = .true.
50			nmflag = .false.
51			step = 0
52
53
54			dynamics: do
55			if (.not.not_done) exit dynamics
56
57			step = step + 1
58
59
60			if (checktemp.and.(mod(step,check_temp_steps).eq.0)) then
61			call calc_temp(verlet_status%temp)
62			if (dabs(verlet_status%temp-target_temp).gt.therm_variance) then
63			call scale_velocities(target_temp)
64			end if
65			end if
66
67			call movea()
68
69
70			call check(update_nlist)
71
72
73			if (do_pot) then
74			call calc_forces(update_nlist,nmflag,pe = verlet_status%pot_e)
75			else
76			call calc_forces(update_nlist,nmflag)
77			endif
78
79			call moveb()
80
81
82			call evolve_time(step,verlet_status,do_pot,not_done)
83
84
85			#ifdef MPI
86			call mpi_barrier(mpi_comm_world,mpi_err)
87			#endif
88			end do dynamics
89			end subroutine velocity_verlet_dynamics
90
91			! *******************************************************************
92			! FICHE F.4
93			! TWO ROUTINES THAT TOGETHER IMPLEMENT VELOCITY VERLET METHOD.
94			!
95			! REFERENCE:
96			!
97			! SWOPE ET AL., J. CHEM. PHYS. 76, 637, 1982.
98			!
99			! ROUTINES SUPPLIED:
100			!
101			! SUBROUTINE MOVEA ( DT, M )
102			! MOVES POSITIONS AND PARTIALLY UPDATES VELOCITIES.
103			! SUBROUTINE MOVEB ( DT, M, K )
104			! COMPLETES VELOCITY MOVE AND CALCULATES KINETIC ENERGY.
105			!
106			! PRINCIPAL VARIABLES:
107			!
108			! INTEGER N NUMBER OF MOLECULES
109			! REAL DT TIMESTEP
110			! REAL M ATOMIC MASS
111			! REAL K KINETIC ENERGY
112			! REAL RX(N),RY(N),RZ(N) POSITIONS
113			! REAL VX(N),VY(N),VZ(N) VELOCITIES
114			! REAL FX(N),FY(N),FZ(N) FORCES
115			!
116			! USAGE:
117			!
118			! AT THE START OF A TIMESTEP, MOVEA IS CALLED TO ADVANCE THE
119			! POSITIONS AND 'HALF-ADVANCE' THE VELOCITIES. THEN THE FORCE
120			! ROUTINE IS CALLED, AND THIS IS FOLLOWED BY MOVEB WHICH
121			! COMPLETES THE ADVANCEMENT OF VELOCITIES.
122			! *******************************************************************
123
124
125
126			subroutine movea( )
127
128			! *******************************************************************
129			! FIRST PART OF VELOCITY VERLET ALGORITHM
130			!
131			! USAGE:
132			!
133			! THE FIRST PART OF THE ALGORITHM IS A TAYLOR SERIES WHICH
134			! ADVANCES POSITIONS FROM T TO T + DT AND VELOCITIES FROM
135			! T TO T + DT/2. AFTER THIS, THE FORCE ROUTINE IS CALLED.
136			! *******************************************************************
137
138
139
140
141			integer :: i, dim
142			!
143			! units for time are femtosec (10^-15 sec)
144			! units for distance are angstroms (10^-10 m)
145			! units for velocity are angstroms femtosec^-1
146			! units for mass are a.m.u. (1.661 * 10^-27 kg)
147			! units for force are kcal mol^-1 angstrom^-1
148			!
149			! converter will put the final terms into angstroms.
150			! or angstrom/femtosecond.
151
152			! converter = 1.0d0/2.390664d3
153
154			! *******************************************************************
155
156
157
158			do i = 1,nlocal
159			do dim = 1, 3
160
161			q(dim,i) = q(dim,i) + dtv(dim,i) + (dtsq2f(dim,i)/mass(i))*KCALMOL_TO_AMUA2FS2
162
163			v(dim,i) = v(dim,i) + (dt2f(dim,i)/mass(i))KCALMOL_TO_AMUA2FS2
164
165			end do
166			end do
167
168			end subroutine movea
169
170
171
172			subroutine moveb()
173			#ifdef MPI
174			use mpi_module
175			#endif
176			! *******************************************************************
177			! SECOND PART OF VELOCITY VERLET ALGORITHM
178			!
179			! USAGE:
180			!
181			! THE SECOND PART OF THE ALGORITHM ADVANCES VELOCITIES FROM
182			! T + DT/2 TO T + DT. THIS ASSUMES THAT FORCES HAVE BEEN
183			! COMPUTED IN THE FORCE ROUTINE AND STORED IN FX, FY, FZ.
184			! *******************************************************************
185
186
187
188			integer :: i, dim
189			!
190			! units for time are femtosec (10^-15 sec)
191			! units for distance are angstroms (10^-10 m)
192			! units for velocity are angstroms femtosec^-1
193			! units for mass are a.m.u. (1.661 * 10^-27 kg)
194			! units for force are kcal mol^-1 angstrom^-1
195			!
196			! converter will put the final terms into kcal/mol
197			! or angstrom/femtosecond.
198
199			! converter = 1.0d0/2.390664d3
200
201			! *******************************************************************
202
203
204
205
206			do i = 1, nlocal
207			do dim = 1, 3
208			v(dim,i) = v(dim,i) + (dt2f(dim,i)/mass(i))KCALMOL_TO_AMUA2FS2
209			end do
210			end do
211
212			end subroutine moveb
213
214			!! modified to be used in a mpi simulation.....
215
216
217
218
219			end module dynamics_velocity_verlet