3 c*********************************************************************72
5 cc MAIN is the main program for MD_OPENMP.
9 c The program implements a simple molecular dynamics simulation.
11 c The program uses Open MP directives to allow parallel computation.
13 c The velocity Verlet time integration scheme is used.
15 c The particles interact with a central pair potential.
17 c Output of the program is saved in the TNG format, which is why this
18 c code is included in the TNG API release.
22 c This code is distributed under the GNU LGPL license.
30 c Original FORTRAN90 version by Bill Magro.
31 c FORTRAN77 version by John Burkardt.
32 c TNG trajectory output by Magnus Lundborg.
45 parameter ( np
= 250 )
47 parameter ( step_num
= 1000 )
49 double precision acc
(nd
,np
)
50 double precision box
(nd
)
51 double precision box_shape
(9)
53 parameter ( dt
= 0.0001D
+00 )
55 double precision force
(nd
,np
)
58 double precision kinetic
60 parameter ( mass
= 1.0D
+00 )
61 double precision pos
(nd
,np
)
62 double precision potential
67 integer step_print_index
68 integer step_print_num
72 double precision vel
(nd
,np
)
73 double precision wtime
76 c Cray pointers are not standard fortran 77, but must be used to allocate
79 pointer
(traj_p
, traj
)
80 pointer
(molecule_p
, molecule
)
81 pointer
(chain_p
, chain
)
82 pointer
(residue_p
, residue
)
83 pointer
(atom_p
, atom
)
91 c The TNG functions expect 64 bit integers
93 integer*8 n_frames_per_frame_set
94 integer*8 frames_saved_cnt
95 integer*8 tng_n_particles
98 c These constants are also defined in tng_io.h, but need to
99 c set in fortran as well. This can be copied to any fortran
100 c source code using the tng_io library.
102 integer*8 TNG_UNCOMPRESSED
103 parameter ( TNG_UNCOMPRESSED
= 0)
104 integer TNG_NON_TRAJECTORY_BLOCK
105 parameter ( TNG_NON_TRAJECTORY_BLOCK
= 0)
106 integer TNG_TRAJECTORY_BLOCK
107 parameter ( TNG_TRAJECTORY_BLOCK
= 1)
108 integer*8 TNG_GENERAL_INFO
109 parameter ( TNG_GENERAL_INFO
= 0 )
110 integer*8 TNG_MOLECULES
111 parameter ( TNG_MOLECULES
= 1 )
112 integer*8 TNG_TRAJECTORY_FRAME_SET
113 parameter ( TNG_TRAJECTORY_FRAME_SET
= 2 )
114 integer*8 TNG_PARTICLE_MAPPING
115 parameter ( TNG_PARTICLE_MAPPING
= 3 )
116 integer*8 TNG_TRAJ_BOX_SHAPE
117 parameter ( TNG_TRAJ_BOX_SHAPE
= 10000 )
118 integer*8 TNG_TRAJ_POSITIONS
119 parameter ( TNG_TRAJ_POSITIONS
= 10001 )
120 integer*8 TNG_TRAJ_VELOCITIES
121 parameter ( TNG_TRAJ_VELOCITIES
= 10002 )
122 integer*8 TNG_TRAJ_FORCES
123 parameter ( TNG_TRAJ_FORCES
= 10003 )
124 integer TNG_SKIP_HASH
125 parameter ( TNG_SKIP_HASH
= 0 )
127 parameter ( TNG_USE_HASH
= 1 )
128 integer TNG_CHAR_DATA
129 parameter ( TNG_CHAR_DATA
= 0 )
131 parameter ( TNG_INT_DATA
= 1 )
132 integer TNG_FLOAT_DATA
133 parameter ( TNG_FLOAT_DATA
= 2 )
134 integer TNG_DOUBLE_DATA
135 parameter ( TNG_DOUBLE_DATA
= 3 )
139 write ( *, '(a)' ) ' '
140 write ( *, '(a)' ) 'MD_OPENMP'
141 write ( *, '(a)' ) ' FORTRAN77/OpenMP version'
142 write ( *, '(a)' ) ' '
143 write ( *, '(a)' ) ' A molecular dynamics program.'
144 write ( *, '(a)' ) ' '
145 write ( *, '(a,i8)' )
146 & ' NP, the number of particles in the simulation is ', np
147 write ( *, '(a,i8)' )
148 & ' STEP_NUM, the number of time steps, is ', step_num
149 write ( *, '(a,g14.6)' )
150 & ' DT, the size of each time step, is ', dt
151 write ( *, '(a)' ) ' '
152 write ( *, '(a,i8)' )
153 & ' The number of processors = ', omp_get_num_procs
( )
154 write ( *, '(a,i8)' )
155 & ' The number of threads = ', omp_get_max_threads
( )
157 write ( *, '(a)' ) ' '
158 write ( *, '(a)' ) ' Initializing trajectory storage.'
159 call tng_trajectory_init
(traj_p
)
162 c N.B. The TNG output file should be modified according to needs
164 call tng_output_file_set
(traj
, TNG_EXAMPLE_FILES_DIR
"tng_md_out_f77.tng")
166 write ( *, '(a)' ) ' Creating molecules in trajectory.'
168 call tng_molecule_add
(traj
, "water", molecule_p
)
169 call tng_molecule_chain_add
(traj
, molecule
, "W", chain_p
)
170 call tng_chain_residue_add
(traj
, chain
, "WAT", residue_p
)
171 call tng_residue_atom_add
(traj
, residue
, "O", "O", atom_p
)
172 call tng_molecule_cnt_set
(traj
, molecule
, tng_n_particles
)
174 c Set the dimensions of the box.
181 box_shape
(i*nd
+ i
) = box
(i
)
185 c Add the box shape data block
187 call tng_data_block_add
(traj
, TNG_TRAJ_BOX_SHAPE
, "BOX SHAPE",
188 & TNG_DOUBLE_DATA
, TNG_NON_TRAJECTORY_BLOCK
, int
(1, 8),
189 & int
(9, 8), int
(1, 8), TNG_UNCOMPRESSED
, box_shape
)
192 c Write the file headers
194 call tng_file_headers_write
(traj
, TNG_USE_HASH
)
197 c Set initial positions, velocities, and accelerations.
200 & ' Initializing positions, velocities, and accelerations.'
202 call initialize
( np
, nd
, box
, seed
, pos
, vel
, acc
)
204 c Compute the forces and energies.
206 write ( *, '(a)' ) ' '
207 write ( *, '(a)' ) ' Computing initial forces and energies.'
209 call compute
( np
, nd
, pos
, vel
, mass
, force
, potential
,
212 e0
= potential
+ kinetic
227 c This is the main time stepping loop.
229 write ( *, '(a,i4,a)' ) ' Every', step_save
,
230 & ' steps particle positions, velocities and forces are'
231 write ( *, '(a)' ) ' saved to a TNG trajectory file.'
234 & ' At each step, we report the potential and kinetic energies.'
236 & ' The sum of these energies should be a constant.'
238 & ' As an accuracy check, we also print the relative error'
239 write ( *, '(a)' ) ' in the total energy.'
240 write ( *, '(a)' ) ' '
242 & ' Step Potential Kinetic (P+K-E0)/E0'
244 & ' Energy P Energy K ' //
245 & 'Relative Energy Error'
246 write ( *, '(a)' ) ' '
249 write ( *, '(2x,i8,2x,g14.6,2x,g14.6,2x,g14.6)' )
250 & step
, potential
, kinetic
, ( potential
+ kinetic
- e0
) / e0
251 step_print_index
= step_print_index
+ 1
252 step_print
= ( step_print_index
* step_num
) / step_print_num
255 c Create a frame set for writing data
257 call tng_num_frames_per_frame_set_get
(traj
,
258 & n_frames_per_frame_set
)
259 call tng_frame_set_new
(traj
, int
(0, 8), n_frames_per_frame_set
)
262 c Add empty data blocks.
264 call tng_particle_data_block_add
(traj
, TNG_TRAJ_POSITIONS
,
265 & "POSITIONS", TNG_DOUBLE_DATA
, TNG_TRAJECTORY_BLOCK
,
266 & n_frames_per_frame_set
, int
(3, 8), int
(1, 8), int
(0, 8),
267 & tng_n_particles
, TNG_UNCOMPRESSED
, %VAL
(int
(0, 8)))
269 call tng_particle_data_block_add
(traj
, TNG_TRAJ_VELOCITIES
,
270 & "VELOCITIES", TNG_DOUBLE_DATA
, TNG_TRAJECTORY_BLOCK
,
271 & n_frames_per_frame_set
, int
(3, 8), int
(1, 8), int
(0, 8),
272 & tng_n_particles
, TNG_UNCOMPRESSED
, %VAL
(int
(0, 8)))
274 call tng_particle_data_block_add
(traj
, TNG_TRAJ_FORCES
,
275 & "FORCES", TNG_DOUBLE_DATA
, TNG_TRAJECTORY_BLOCK
,
276 & n_frames_per_frame_set
, int
(3, 8), int
(1, 8), int
(0, 8),
277 & tng_n_particles
, TNG_UNCOMPRESSED
, %VAL
(int
(0, 8)))
280 c The potential energy data block is saved sparsely.
282 call tng_data_block_add
(traj
, int
(10101, 8),
283 & "POTENTIAL ENERGY", TNG_DOUBLE_DATA
, TNG_TRAJECTORY_BLOCK
,
284 & n_frames_per_frame_set
, int
(1, 8), sparse_save
,
285 & TNG_UNCOMPRESSED
, %VAL
(int
(0, 8)))
289 c Write the frame set to disk
291 call tng_frame_set_write
(traj
, TNG_USE_HASH
)
293 wtime
= omp_get_wtime
( )
295 do step
= 1, step_num
297 call compute
( np
, nd
, pos
, vel
, mass
, force
, potential
,
300 if ( step
.eq
. step_print
) then
302 write ( *, '(2x,i8,2x,g14.6,2x,g14.6,2x,g14.6)' )
303 & step
, potential
, kinetic
, ( potential
+ kinetic
- e0
) / e0
305 step_print_index
= step_print_index
+ 1
306 step_print
= ( step_print_index
* step_num
) / step_print_num
311 c Output to TNG file at regular intervals, specified by step_save
313 if ( step_save
.EQ
. 0 .OR
. mod
(step
, step_save
) .EQ
. 0 ) then
314 call tng_frame_particle_data_write
(traj
, frames_saved_cnt
,
315 & TNG_TRAJ_POSITIONS
, int
(0, 8), tng_n_particles
, pos
,
317 call tng_frame_particle_data_write
(traj
, frames_saved_cnt
,
318 & TNG_TRAJ_VELOCITIES
, int
(0, 8), tng_n_particles
, vel
,
320 call tng_frame_particle_data_write
(traj
, frames_saved_cnt
,
321 & TNG_TRAJ_FORCES
, int
(0, 8), tng_n_particles
, force
,
323 frames_saved_cnt
= frames_saved_cnt
+ 1
325 if (mod
(step
, step_save
* sparse_save
) .EQ
. 0) then
326 call tng_frame_data_write
(traj
, frames_saved_cnt
,
327 & int
(10101, 8), potential
, TNG_USE_HASH
)
332 call update
( np
, nd
, pos
, vel
, force
, acc
, mass
, dt
)
336 wtime
= omp_get_wtime
( ) - wtime
338 write ( *, '(a)' ) ' '
340 & ' Elapsed time for main computation:'
341 write ( *, '(2x,g14.6,a)' ) wtime
, ' seconds'
345 call tng_trajectory_destroy
(traj_p
)
347 write ( *, '(a)' ) ' '
348 write ( *, '(a)' ) 'MD_OPENMP'
349 write ( *, '(a)' ) ' Normal end of execution.'
351 write ( *, '(a)' ) ' '
356 subroutine compute
( np
, nd
, pos
, vel
, mass
, f
, pot
, kin
)
358 c*********************************************************************72
360 cc COMPUTE computes the forces and energies.
364 c The computation of forces and energies is fully parallel.
368 c This code is distributed under the GNU LGPL license.
376 c Original FORTRAN90 version by Bill Magro.
377 c FORTRAN77 version by John Burkardt.
381 c Input, integer NP, the number of particles.
383 c Input, integer ND, the number of spatial dimensions.
385 c Input, double precision POS(ND,NP), the position of each particle.
387 c Input, double precision VEL(ND,NP), the velocity of each particle.
389 c Input, double precision MASS, the mass of each particle.
399 double precision f
(nd
,np
)
404 double precision mass
406 parameter ( PI2
= 3.141592653589793D
+00 / 2.0D
+00 )
407 double precision pos
(nd
,np
)
409 double precision rij
(nd
)
411 double precision vel
(nd
,np
)
417 c$omp& shared ( f, nd, np, pos, vel )
418 c$omp& private ( d, d2, i, j, k, rij )
420 c$omp do reduction ( + : pot, kin )
423 c Compute the potential energy and forces.
433 call dist
( nd
, pos
(1,i
), pos
(1,j
), rij
, d
)
435 c Attribute half of the potential energy to particle J.
439 pot
= pot
+ 0.5D
+00 * ( sin
( d2
) )**2
442 f
(k
,i
) = f
(k
,i
) - rij
(k
) * sin
( 2.0D
+00 * d2
) / d
449 c Compute the kinetic energy.
452 kin
= kin
+ vel
(k
,i
)**2
460 kin
= kin
* 0.5D
+00 * mass
464 subroutine dist
( nd
, r1
, r2
, dr
, d
)
466 c*********************************************************************72
468 cc DIST computes the displacement (and its norm) between two particles.
472 c This code is distributed under the GNU LGPL license.
480 c Original FORTRAN90 version by Bill Magro.
481 c FORTRAN77 version by John Burkardt.
485 c Input, integer ND, the number of spatial dimensions.
487 c Input, double precision R1(ND), R2(ND), the positions of the particles.
489 c Output, double precision DR(ND), the displacement vector.
491 c Output, double precision D, the Euclidean norm of the displacement.
498 double precision dr
(nd
)
500 double precision r1
(nd
)
501 double precision r2
(nd
)
504 dr
(i
) = r1
(i
) - r2
(i
)
515 subroutine initialize
( np
, nd
, box
, seed
, pos
, vel
, acc
)
517 c*********************************************************************72
519 cc INITIALIZE initializes the positions, velocities, and accelerations.
523 c This code is distributed under the GNU LGPL license.
531 c Original FORTRAN90 version by Bill Magro.
532 c FORTRAN77 version by John Burkardt.
536 c Input, integer NP, the number of particles.
538 c Input, integer ND, the number of spatial dimensions.
540 c Input, double precision BOX(ND), specifies the maximum position
541 c of particles in each dimension.
543 c Input/output, integer SEED, a seed for the random number generator.
545 c Output, double precision POS(ND,NP), the position of each particle.
547 c Output, double precision VEL(ND,NP), the velocity of each particle.
549 c Output, double precision ACC(ND,NP), the acceleration of each particle.
556 double precision acc
(nd
,np
)
557 double precision box
(nd
)
560 double precision pos
(nd
,np
)
561 double precision r8_uniform_01
563 double precision vel
(nd
,np
)
565 c Give the particles random positions within the box.
569 pos
(i
,j
) = r8_uniform_01
( seed
)
574 c$omp& shared ( acc, box, nd, np, pos, vel )
575 c$omp& private ( i, j )
580 pos
(i
,j
) = box
(i
) * pos
(i
,j
)
591 function r8_uniform_01
( seed
)
593 c*********************************************************************72
595 cc R8_UNIFORM_01 returns a unit pseudorandom R8.
599 c This routine implements the recursion
601 c seed = 16807 * seed mod ( 2**31 - 1 )
602 c r8_uniform_01 = seed / ( 2**31 - 1 )
604 c The integer arithmetic never requires more than 32 bits,
605 c including a sign bit.
607 c If the initial seed is 12345, then the first three computations are
609 c Input Output R8_UNIFORM_01
612 c 12345 207482415 0.096616
613 c 207482415 1790989824 0.833995
614 c 1790989824 2035175616 0.947702
618 c This code is distributed under the GNU LGPL license.
630 c Paul Bratley, Bennett Fox, Linus Schrage,
631 c A Guide to Simulation,
632 c Springer Verlag, pages 201-202, 1983.
635 c Random Number Generation,
636 c in Handbook of Simulation,
637 c edited by Jerry Banks,
638 c Wiley Interscience, page 95, 1998.
642 c Implementation and Relative Efficiency of Quasirandom
643 c Sequence Generators,
644 c ACM Transactions on Mathematical Software,
645 c Volume 12, Number 4, pages 362-376, 1986.
647 c Peter Lewis, Allen Goodman, James Miller,
648 c A Pseudo-Random Number Generator for the System/360,
649 c IBM Systems Journal,
650 c Volume 8, pages 136-143, 1969.
654 c Input/output, integer SEED, the "seed" value, which should NOT be 0.
655 c On output, SEED has been updated.
657 c Output, double precision R8_UNIFORM_01, a new pseudorandom variate,
658 c strictly between 0 and 1.
662 double precision r8_uniform_01
666 if ( seed
.eq
. 0 ) then
667 write ( *, '(a)' ) ' '
668 write ( *, '(a)' ) 'R8_UNIFORM_01 - Fatal error!'
669 write ( *, '(a)' ) ' Input value of SEED = 0.'
675 seed
= 16807 * ( seed
- k
* 127773 ) - k
* 2836
677 if ( seed
.lt
. 0 ) then
678 seed
= seed
+ 2147483647
681 c Although SEED can be represented exactly as a 32 bit integer,
682 c it generally cannot be represented exactly as a 32 bit real number!
684 r8_uniform_01
= dble
( seed
) * 4.656612875D
-10
688 subroutine timestamp
( )
690 c*********************************************************************72
692 cc TIMESTAMP prints out the current YMDHMS date as a timestamp.
696 c This code is distributed under the GNU LGPL license.
712 character * ( 8 ) ampm
714 character * ( 8 ) date
718 character * ( 9 ) month
(12)
721 character * ( 10 ) time
727 & 'January ', 'February ', 'March ', 'April ',
728 & 'May ', 'June ', 'July ', 'August ',
729 & 'September', 'October ', 'November ', 'December ' /
731 call date_and_time
( date, time
)
733 read ( date, '(i4,i2,i2)' ) y
, m
, d
734 read ( time
, '(i2,i2,i2,1x,i3)' ) h
, n
, s
, mm
736 if ( h
.lt
. 12 ) then
738 else if ( h
.eq
. 12 ) then
739 if ( n
.eq
. 0 .and
. s
.eq
. 0 ) then
746 if ( h
.lt
. 12 ) then
748 else if ( h
.eq
. 12 ) then
749 if ( n
.eq
. 0 .and
. s
.eq
. 0 ) then
758 & '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' )
759 & d
, month
(m
), y
, h
, ':', n
, ':', s
, '.', mm
, ampm
763 subroutine update
( np
, nd
, pos
, vel
, f
, acc
, mass
, dt
)
765 c*********************************************************************72
767 cc UPDATE performs the time integration, using a velocity Verlet algorithm.
771 c The time integration is fully parallel.
775 c This code is distributed under the GNU LGPL license.
783 c Original FORTRAN90 version by Bill Magro.
784 c FORTRAN77 version by John Burkardt.
788 c Input, integer NP, the number of particles.
790 c Input, integer ND, the number of spatial dimensions.
792 c Input/output, double precision POS(ND,NP), the position of each particle.
794 c Input/output, double precision VEL(ND,NP), the velocity of each particle.
796 c Input, double precision MASS, the mass of each particle.
798 c Input/output, double precision ACC(ND,NP), the acceleration of each
806 double precision acc
(nd
,np
)
808 double precision f
(nd
,np
)
811 double precision mass
812 double precision pos
(nd
,np
)
813 double precision rmass
814 double precision vel
(nd
,np
)
816 rmass
= 1.0D
+00 / mass
819 c$omp& shared ( acc, dt, f, nd, np, pos, rmass, vel )
820 c$omp& private ( i, j )
827 & + vel
(i
,j
) * dt
+ 0.5D
+00 * acc
(i
,j
) * dt
* dt
830 & + 0.5D
+00 * dt
* ( f
(i
,j
) * rmass
+ acc
(i
,j
) )
832 acc
(i
,j
) = f
(i
,j
) * rmass