1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     #################################################################
9c     ##                                                             ##
10c     ##  program dynamic  --  run molecular or stochastic dynamics  ##
11c     ##                                                             ##
12c     #################################################################
13c
14c
15c     "dynamic" computes a molecular or stochastic dynamics trajectory
16c     in one of the standard statistical mechanical ensembles and using
17c     any of several possible integration methods
18c
19c
20      program dynamic
21      use atoms
22      use bath
23      use bndstr
24      use bound
25      use inform
26      use iounit
27      use keys
28      use mdstuf
29      use potent
30      use stodyn
31      use usage
32      implicit none
33      integer i,istep,nstep
34      integer mode,next
35      real*8 dt,dtsave
36      logical exist
37      character*20 keyword
38      character*240 record
39      character*240 string
40c
41c
42c     set up the structure and molecular mechanics calculation
43c
44      call initial
45      call getxyz
46      call mechanic
47c
48c     initialize the temperature, pressure and coupling baths
49c
50      kelvin = 0.0d0
51      atmsph = 0.0d0
52      isothermal = .false.
53      isobaric = .false.
54c
55c     check for keywords containing any altered parameters
56c
57      integrate = 'BEEMAN'
58      do i = 1, nkey
59         next = 1
60         record = keyline(i)
61         call gettext (record,keyword,next)
62         call upcase (keyword)
63         string = record(next:240)
64         if (keyword(1:11) .eq. 'INTEGRATOR ') then
65            call getword (record,integrate,next)
66            call upcase (integrate)
67         end if
68      end do
69c
70c     initialize the simulation length as number of time steps
71c
72      nstep = -1
73      call nextarg (string,exist)
74      if (exist)  read (string,*,err=10,end=10)  nstep
75   10 continue
76      do while (nstep .lt. 0)
77         write (iout,20)
78   20    format (/,' Enter the Number of Dynamics Steps to be',
79     &              ' Taken :  ',$)
80         read (input,30,err=40)  nstep
81   30    format (i10)
82         if (nstep .lt. 0)  nstep = 0
83   40    continue
84      end do
85c
86c     get the length of the dynamics time step in picoseconds
87c
88      dt = -1.0d0
89      call nextarg (string,exist)
90      if (exist)  read (string,*,err=50,end=50)  dt
91   50 continue
92      do while (dt .lt. 0.0d0)
93         write (iout,60)
94   60    format (/,' Enter the Time Step Length in Femtoseconds',
95     &              ' [1.0] :  ',$)
96         read (input,70,err=80)  dt
97   70    format (f20.0)
98         if (dt .le. 0.0d0)  dt = 1.0d0
99   80    continue
100      end do
101      dt = 0.001d0 * dt
102c
103c     enforce bounds on thermostat and barostat coupling times
104c
105      tautemp = max(tautemp,dt)
106      taupres = max(taupres,dt)
107c
108c     set the time between trajectory snapshot coordinate saves
109c
110      dtsave = -1.0d0
111      call nextarg (string,exist)
112      if (exist)  read (string,*,err=90,end=90)  dtsave
113   90 continue
114      do while (dtsave .lt. 0.0d0)
115         write (iout,100)
116  100    format (/,' Enter Time between Saves in Picoseconds',
117     &              ' [0.1] :  ',$)
118         read (input,110,err=120)  dtsave
119  110    format (f20.0)
120         if (dtsave .le. 0.0d0)  dtsave = 0.1d0
121  120    continue
122      end do
123      iwrite = nint(dtsave/dt)
124c
125c     get choice of statistical ensemble for periodic system
126c
127      if (use_bounds) then
128         mode = -1
129         call nextarg (string,exist)
130         if (exist)  read (string,*,err=130,end=130)  mode
131  130    continue
132         do while (mode.lt.1 .or. mode.gt.4)
133            write (iout,140)
134  140       format (/,' Available Statistical Mechanical Ensembles :',
135     &              //,4x,'(1) Microcanonical (NVE)',
136     &              /,4x,'(2) Canonical (NVT)',
137     &              /,4x,'(3) Isoenthalpic-Isobaric (NPH)',
138     &              /,4x,'(4) Isothermal-Isobaric (NPT)',
139     &              //,' Enter the Number of the Desired Choice',
140     &                 ' [1] :  ',$)
141            read (input,150,err=160)  mode
142  150       format (i10)
143            if (mode .le. 0)  mode = 1
144  160       continue
145         end do
146         if (integrate.eq.'BUSSI' .or. integrate.eq.'NOSE-HOOVER'
147     &                .or. integrate.eq.'GHMC') then
148            if (mode .ne. 4) then
149               mode = 4
150               write (iout,170)
151  170          format (/,' Switching to NPT Ensemble as Required',
152     &                    ' by Chosen Integrator')
153            end if
154         end if
155         if (mode.eq.2 .or. mode.eq.4) then
156            isothermal = .true.
157            kelvin = -1.0d0
158            call nextarg (string,exist)
159            if (exist)  read (string,*,err=180,end=180)  kelvin
160  180       continue
161            do while (kelvin .le. 0.0d0)
162               write (iout,190)
163  190          format (/,' Enter the Desired Temperature in Degrees',
164     &                    ' K [298] :  ',$)
165               read (input,200,err=210)  kelvin
166  200          format (f20.0)
167               if (kelvin .le. 0.0d0)  kelvin = 298.0d0
168  210          continue
169            end do
170         end if
171         if (mode.eq.3 .or. mode.eq.4) then
172            isobaric = .true.
173            atmsph = -1.0d0
174            call nextarg (string,exist)
175            if (exist)  read (string,*,err=220,end=220)  atmsph
176  220       continue
177            do while (atmsph .eq. -1.0d0)
178               write (iout,230)
179  230          format (/,' Enter the Desired Pressure in Atm',
180     &                    ' [1.0] :  ',$)
181               read (input,240,err=250)  atmsph
182  240          format (f20.0)
183               if (atmsph .eq. -1.0d0)  atmsph = 1.0d0
184  250          continue
185            end do
186         end if
187      end if
188c
189c     use constant energy or temperature for nonperiodic system
190c
191      if (.not. use_bounds) then
192         mode = -1
193         call nextarg (string,exist)
194         if (exist)  read (string,*,err=260,end=260)  mode
195  260    continue
196         do while (mode.lt.1 .or. mode.gt.2)
197            write (iout,270)
198  270       format (/,' Available Simulation Control Modes :',
199     &              //,4x,'(1) Constant Total Energy Value (E)',
200     &              /,4x,'(2) Constant Temperature via Thermostat (T)',
201     &              //,' Enter the Number of the Desired Choice',
202     &                 ' [1] :  ',$)
203            read (input,280,err=290)  mode
204  280       format (i10)
205            if (mode .le. 0)  mode = 1
206  290       continue
207         end do
208         if (mode .eq. 2) then
209            isothermal = .true.
210            kelvin = -1.0d0
211            call nextarg (string,exist)
212            if (exist)  read (string,*,err=300,end=300)  kelvin
213  300       continue
214            do while (kelvin .le. 0.0d0)
215               write (iout,310)
216  310          format (/,' Enter the Desired Temperature in Degrees',
217     &                    ' K [298] :  ',$)
218               read (input,320,err=330)  kelvin
219  320          format (f20.0)
220               if (kelvin .le. 0.0d0)  kelvin = 298.0d0
221  330          continue
222            end do
223         end if
224      end if
225c
226c     perform the setup functions needed to run dynamics
227c
228      call mdinit
229c
230c     print out a header line for the dynamics computation
231c
232      if (integrate .eq. 'VERLET') then
233         write (iout,340)
234  340    format (/,' Molecular Dynamics Trajectory via',
235     &              ' Velocity Verlet Algorithm')
236      else if (integrate .eq. 'STOCHASTIC') then
237         write (iout,350)
238  350    format (/,' Stochastic Dynamics Trajectory via',
239     &              ' Velocity Verlet Algorithm')
240      else if (integrate .eq. 'BAOAB') then
241         write (iout,360)
242  360    format (/,' Constrained Stochastic Dynamics Trajectory',
243     &              ' via BAOAB Algorithm')
244      else if (integrate .eq. 'BUSSI') then
245         write (iout,370)
246  370    format (/,' Molecular Dynamics Trajectory via',
247     &              ' Bussi-Parrinello NPT Algorithm')
248      else if (integrate .eq. 'NOSE-HOOVER') then
249         write (iout,380)
250  380    format (/,' Molecular Dynamics Trajectory via',
251     &              ' Nose-Hoover NPT Algorithm')
252      else if (integrate .eq. 'GHMC') then
253         write (iout,390)
254  390    format (/,' Stochastic Dynamics Trajectory via',
255     &              ' Generalized Hybrid Monte Carlo')
256      else if (integrate .eq. 'RIGIDBODY') then
257         write (iout,400)
258  400    format (/,' Molecular Dynamics Trajectory via',
259     &              ' Rigid Body Algorithm')
260      else if (integrate .eq. 'RESPA') then
261         write (iout,410)
262  410    format (/,' Molecular Dynamics Trajectory via',
263     &              ' r-RESPA MTS Algorithm')
264      else
265         write (iout,420)
266  420    format (/,' Molecular Dynamics Trajectory via',
267     &              ' Modified Beeman Algorithm')
268      end if
269      flush (iout)
270c
271c     integrate equations of motion to take a time step
272c
273      do istep = 1, nstep
274         if (integrate .eq. 'VERLET') then
275            call verlet (istep,dt)
276         else if (integrate .eq. 'STOCHASTIC') then
277            call sdstep (istep,dt)
278         else if (integrate .eq. 'BAOAB') then
279            call baoab (istep,dt)
280         else if (integrate .eq. 'BUSSI') then
281            call bussi (istep,dt)
282         else if (integrate .eq. 'NOSE-HOOVER') then
283            call nose (istep,dt)
284         else if (integrate .eq. 'GHMC') then
285            call ghmcstep (istep,dt)
286         else if (integrate .eq. 'RIGIDBODY') then
287            call rgdstep (istep,dt)
288         else if (integrate .eq. 'RESPA') then
289            call respa (istep,dt)
290         else
291            call beeman (istep,dt)
292         end if
293      end do
294c
295c     perform any final tasks before program exit
296c
297      call final
298      end
299