wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / source / unres / src_MD-M / MD.F
1       subroutine MD
2 c------------------------------------------------
3 c  The driver for molecular dynamics subroutines
4 c------------------------------------------------
5       implicit real*8 (a-h,o-z)
6       include 'DIMENSIONS'
7 #ifdef MPI
8       include "mpif.h"
9 #endif
10       include 'COMMON.CONTROL'
11       include 'COMMON.VAR'
12       include 'COMMON.MD'
13 #ifndef LANG0
14       include 'COMMON.LANGEVIN'
15 #else
16       include 'COMMON.LANGEVIN.lang0'
17 #endif
18       include 'COMMON.CHAIN'
19       include 'COMMON.DERIV'
20       include 'COMMON.GEO'
21       include 'COMMON.LOCAL'
22       include 'COMMON.INTERACT'
23       include 'COMMON.IOUNITS'
24       include 'COMMON.NAMES'
25       include 'COMMON.TIME1'
26       include 'COMMON.HAIRPIN'
27       double precision cm(3),L(3),vcm(3)
28 #ifdef VOUT
29       double precision v_work(maxres6),v_transf(maxres6)
30 #endif
31       integer ilen,rstcount
32       external ilen
33       character*50 tytul
34       common /gucio/ cm
35       integer itime
36 c
37 #ifdef MPI
38       if (ilen(tmpdir).gt.0)
39      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
40      &        //liczba(:ilen(liczba))//'.rst')
41 #else
42       if (ilen(tmpdir).gt.0)
43      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
44 #endif
45       t_MDsetup=0.0d0
46       t_langsetup=0.0d0
47       t_MD=0.0d0
48       t_enegrad=0.0d0
49       t_sdsetup=0.0d0
50       write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
51 #ifdef MPI
52       tt0=MPI_Wtime()
53 #else
54       tt0 = tcpu()
55 #endif
56 c Determine the inverse of the inertia matrix.
57       call setup_MD_matrices
58 c Initialize MD
59       call init_MD
60 #ifdef MPI
61       t_MDsetup = MPI_Wtime()-tt0
62 #else
63       t_MDsetup = tcpu()-tt0
64 #endif
65       rstcount=0 
66 c   Entering the MD loop       
67 #ifdef MPI
68       tt0 = MPI_Wtime()
69 #else
70       tt0 = tcpu()
71 #endif
72       if (lang.eq.2 .or. lang.eq.3) then
73         call setup_fricmat
74         if (lang.eq.2) then
75           call sd_verlet_p_setup        
76         else
77           call sd_verlet_ciccotti_setup
78         endif
79 #ifndef   LANG0
80         do i=1,dimen
81           do j=1,dimen
82             pfric0_mat(i,j,0)=pfric_mat(i,j)
83             afric0_mat(i,j,0)=afric_mat(i,j)
84             vfric0_mat(i,j,0)=vfric_mat(i,j)
85             prand0_mat(i,j,0)=prand_mat(i,j)
86             vrand0_mat1(i,j,0)=vrand_mat1(i,j)
87             vrand0_mat2(i,j,0)=vrand_mat2(i,j)
88           enddo
89         enddo
90 #endif
91         flag_stoch(0)=.true.
92         do i=1,maxflag_stoch
93           flag_stoch(i)=.false.
94         enddo  
95       else if (lang.eq.1 .or. lang.eq.4) then
96         call setup_fricmat
97       endif
98 #ifdef MPI
99       t_langsetup=MPI_Wtime()-tt0
100       tt0=MPI_Wtime()
101 #else
102       t_langsetup=tcpu()-tt0
103       tt0=tcpu()
104 #endif
105       do itime=1,n_timestep
106         if (ovrtim()) goto 100
107         rstcount=rstcount+1
108         if (lang.gt.0 .and. surfarea .and. 
109      &      mod(itime,reset_fricmat).eq.0) then
110           if (lang.eq.2 .or. lang.eq.3) then
111             call setup_fricmat
112             if (lang.eq.2) then
113               call sd_verlet_p_setup
114             else
115               call sd_verlet_ciccotti_setup
116             endif
117 #ifndef   LANG0
118             do i=1,dimen
119               do j=1,dimen
120                 pfric0_mat(i,j,0)=pfric_mat(i,j)
121                 afric0_mat(i,j,0)=afric_mat(i,j)
122                 vfric0_mat(i,j,0)=vfric_mat(i,j)
123                 prand0_mat(i,j,0)=prand_mat(i,j)
124                 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
125                 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
126               enddo
127             enddo
128 #endif
129             flag_stoch(0)=.true.
130             do i=1,maxflag_stoch
131               flag_stoch(i)=.false.
132             enddo   
133           else if (lang.eq.1 .or. lang.eq.4) then
134             call setup_fricmat
135           endif
136           write (iout,'(a,i10)') 
137      &      "Friction matrix reset based on surface area, itime",itime
138         endif
139         if (reset_vel .and. tbf .and. lang.eq.0 
140      &      .and. mod(itime,count_reset_vel).eq.0) then
141           call random_vel
142           write(iout,'(a,f20.2)') 
143      &     "Velocities reset to random values, time",totT       
144           do i=0,2*nres
145             do j=1,3
146               d_t_old(j,i)=d_t(j,i)
147             enddo
148           enddo
149         endif
150         if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
151           call inertia_tensor  
152           call vcm_vel(vcm)
153           do j=1,3
154              d_t(j,0)=d_t(j,0)-vcm(j)
155           enddo
156           call kinetic(EK)
157           kinetic_T=2.0d0/(dimen*Rb)*EK
158           scalfac=dsqrt(T_bath/kinetic_T)
159           write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT       
160           do i=0,2*nres
161             do j=1,3
162               d_t_old(j,i)=scalfac*d_t(j,i)
163             enddo
164           enddo
165         endif  
166         if (lang.ne.4) then
167           if (RESPA) then
168 c Time-reversible RESPA algorithm 
169 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
170             call RESPA_step(itime)
171           else
172 c Variable time step algorithm.
173             call velverlet_step(itime)
174           endif
175         else
176           call brown_step(itime)
177         endif
178         if (ntwe.ne.0) then
179          if (mod(itime,ntwe).eq.0) then
180            call statout(itime)
181            call returnbox
182           endif
183 #ifdef VOUT
184         do j=1,3
185           v_work(j)=d_t(j,0)
186         enddo
187         ind=3
188         do i=nnt,nct-1
189           do j=1,3
190             ind=ind+1
191             v_work(ind)=d_t(j,i)
192           enddo
193         enddo
194         do i=nnt,nct
195           if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
196             do j=1,3
197               ind=ind+1
198               v_work(ind)=d_t(j,i+nres)
199             enddo
200           endif
201         enddo
202
203         write (66,'(80f10.5)') 
204      &    ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
205         do i=1,ind
206           v_transf(i)=0.0d0
207           do j=1,ind
208             v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
209           enddo
210            v_transf(i)= v_transf(i)*dsqrt(geigen(i))
211         enddo
212         write (67,'(80f10.5)') (v_transf(i),i=1,ind)
213 #endif
214         endif
215         if (mod(itime,ntwx).eq.0) then
216           write (tytul,'("time",f8.2)') totT
217           if(mdpdb) then
218              call hairpin(.true.,nharp,iharp)
219              call secondary2(.true.)
220              call pdbout(potE,tytul,ipdb)
221           else 
222              call cartout(totT)
223           endif
224         endif
225         if (rstcount.eq.1000.or.itime.eq.n_timestep) then
226            open(irest2,file=rest2name,status='unknown')
227            write(irest2,*) totT,EK,potE,totE,t_bath
228            do i=1,2*nres
229             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
230            enddo
231            do i=1,2*nres
232             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
233            enddo
234           close(irest2)
235           rstcount=0
236         endif 
237       enddo
238   100 continue
239 #ifdef MPI
240       t_MD=MPI_Wtime()-tt0
241 #else
242       t_MD=tcpu()-tt0
243 #endif
244       write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') 
245      &  '  Timing  ',
246      & 'MD calculations setup:',t_MDsetup,
247      & 'Energy & gradient evaluation:',t_enegrad,
248      & 'Stochastic MD setup:',t_langsetup,
249      & 'Stochastic MD step setup:',t_sdsetup,
250      & 'MD steps:',t_MD
251       write (iout,'(/28(1h=),a25,27(1h=))') 
252      & '  End of MD calculation  '
253       return
254       end  
255 c-------------------------------------------------------------------------------
256       subroutine brown_step(itime)
257 c------------------------------------------------
258 c  Perform a single Euler integration step of Brownian dynamics
259 c------------------------------------------------
260       implicit real*8 (a-h,o-z)
261       include 'DIMENSIONS'
262 #ifdef MPI
263       include 'mpif.h'
264 #endif
265       include 'COMMON.CONTROL'
266       include 'COMMON.VAR'
267       include 'COMMON.MD'
268 #ifndef LANG0
269       include 'COMMON.LANGEVIN'
270 #else
271       include 'COMMON.LANGEVIN.lang0'
272 #endif
273       include 'COMMON.CHAIN'
274       include 'COMMON.DERIV'
275       include 'COMMON.GEO'
276       include 'COMMON.LOCAL'
277       include 'COMMON.INTERACT'
278       include 'COMMON.IOUNITS'
279       include 'COMMON.NAMES'
280       include 'COMMON.TIME1'
281       double precision zapas(MAXRES6)
282       integer ilen,rstcount
283       external ilen
284       double precision stochforcvec(MAXRES6)
285       double precision Bmat(MAXRES6,MAXRES2),Cmat(maxres2,maxres2),
286      & Cinv(maxres2,maxres2),GBmat(MAXRES6,MAXRES2),
287      & Tmat(MAXRES6,MAXRES2),Pmat(maxres6,maxres6),Td(maxres6),
288      & ppvec(maxres2)
289       common /stochcalc/ stochforcvec
290       common /gucio/ cm
291       integer itime
292       logical lprn /.false./,lprn1 /.false./
293       integer maxiter /5/
294       double precision difftol /1.0d-5/
295       nbond=nct-nnt
296       do i=nnt,nct
297         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) nbond=nbond+1
298       enddo
299 c
300       if (lprn1) then
301         write (iout,*) "Generalized inverse of fricmat"
302         call matout(dimen,dimen,MAXRES6,MAXRES6,fricmat)
303       endif 
304       do i=1,dimen
305         do j=1,nbond
306           Bmat(i,j)=0.0d0
307         enddo
308       enddo
309       ind=3
310       ind1=0
311       do i=nnt,nct-1
312         ind1=ind1+1
313         do j=1,3
314           Bmat(ind+j,ind1)=dC_norm(j,i)
315         enddo
316         ind=ind+3
317       enddo
318       do i=nnt,nct
319         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
320           ind1=ind1+1
321           do j=1,3
322             Bmat(ind+j,ind1)=dC_norm(j,i+nres)
323           enddo
324           ind=ind+3
325         endif
326       enddo
327       if (lprn1) then 
328         write (iout,*) "Matrix Bmat"
329         call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Bmat)
330       endif
331       do i=1,dimen
332         do j=1,nbond
333           GBmat(i,j)=0.0d0
334           do k=1,dimen
335             GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
336           enddo
337         enddo
338       enddo   
339       if (lprn1) then
340         write (iout,*) "Matrix GBmat"
341         call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Gbmat)
342       endif
343       do i=1,nbond
344         do j=1,nbond
345           Cmat(i,j)=0.0d0
346           do k=1,dimen
347             Cmat(i,j)=Cmat(i,j)+Bmat(k,i)*GBmat(k,j)
348           enddo
349         enddo
350       enddo
351       if (lprn1) then
352         write (iout,*) "Matrix Cmat"
353         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
354       endif
355       call matinvert(nbond,MAXRES2,Cmat,Cinv) 
356       if (lprn1) then
357         write (iout,*) "Matrix Cinv"
358         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cinv)
359       endif
360       do i=1,dimen
361         do j=1,nbond
362           Tmat(i,j)=0.0d0
363           do k=1,nbond
364             Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
365           enddo
366         enddo
367       enddo
368       if (lprn1) then
369         write (iout,*) "Matrix Tmat"
370         call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Tmat)
371       endif
372       do i=1,dimen
373         do j=1,dimen
374           if (i.eq.j) then
375             Pmat(i,j)=1.0d0
376           else
377             Pmat(i,j)=0.0d0
378           endif
379           do k=1,nbond
380             Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
381           enddo
382         enddo
383       enddo
384       if (lprn1) then
385         write (iout,*) "Matrix Pmat"
386         call MATOUT(dimen,dimen,MAXRES6,MAXRES6,Pmat)
387       endif
388       do i=1,dimen
389         Td(i)=0.0d0
390         ind=0
391         do k=nnt,nct-1
392           ind=ind+1
393           Td(i)=Td(i)+vbl*Tmat(i,ind)
394         enddo
395         do k=nnt,nct
396           if (itype(k).ne.10 .and. itype(i).ne.ntyp1) then
397             ind=ind+1
398             Td(i)=Td(i)+vbldsc0(1,itype(k))*Tmat(i,ind)
399           endif
400         enddo
401       enddo 
402       if (lprn1) then
403         write (iout,*) "Vector Td"
404         do i=1,dimen
405           write (iout,'(i5,f10.5)') i,Td(i)
406         enddo
407       endif
408       call stochastic_force(stochforcvec)
409       if (lprn) then
410         write (iout,*) "stochforcvec"
411         do i=1,dimen
412           write (iout,*) i,stochforcvec(i)
413         enddo
414       endif
415       do j=1,3
416         zapas(j)=-gcart(j,0)+stochforcvec(j)
417         d_t_work(j)=d_t(j,0)
418         dC_work(j)=dC_old(j,0)
419       enddo
420       ind=3      
421       do i=nnt,nct-1
422         do j=1,3
423           ind=ind+1
424           zapas(ind)=-gcart(j,i)+stochforcvec(ind)
425           dC_work(ind)=dC_old(j,i)
426         enddo
427       enddo
428       do i=nnt,nct
429         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
430           do j=1,3
431             ind=ind+1
432             zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
433             dC_work(ind)=dC_old(j,i+nres)
434           enddo
435         endif
436       enddo
437
438       if (lprn) then
439         write (iout,*) "Initial d_t_work"
440         do i=1,dimen
441           write (iout,*) i,d_t_work(i)
442         enddo
443       endif
444
445       do i=1,dimen
446         d_t_work(i)=0.0d0
447         do j=1,dimen
448           d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
449         enddo
450       enddo
451
452       do i=1,dimen
453         zapas(i)=Td(i)
454         do j=1,dimen
455           zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
456         enddo
457       enddo
458       if (lprn1) then
459         write (iout,*) "Final d_t_work and zapas"
460         do i=1,dimen
461           write (iout,*) i,d_t_work(i),zapas(i)
462         enddo
463       endif
464
465       do j=1,3
466         d_t(j,0)=d_t_work(j)
467         dc(j,0)=zapas(j)
468         dc_work(j)=dc(j,0)
469       enddo
470       ind=3
471       do i=nnt,nct-1
472         do j=1,3
473           d_t(j,i)=d_t_work(i)
474           dc(j,i)=zapas(ind+j)
475           dc_work(ind+j)=dc(j,i)
476         enddo
477         ind=ind+3
478       enddo
479       do i=nnt,nct
480         do j=1,3
481           d_t(j,i+nres)=d_t_work(ind+j)
482           dc(j,i+nres)=zapas(ind+j)
483           dc_work(ind+j)=dc(j,i+nres)
484         enddo
485         ind=ind+3
486       enddo
487       if (lprn) then
488         call chainbuild_cart
489         write (iout,*) "Before correction for rotational lengthening"
490         write (iout,*) "New coordinates",
491      &  " and differences between actual and standard bond lengths"
492         ind=0
493         do i=nnt,nct-1
494           ind=ind+1
495           xx=vbld(i+1)-vbl
496           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
497      &        i,(dC(j,i),j=1,3),xx
498         enddo
499         do i=nnt,nct
500           if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
501             ind=ind+1
502             xx=vbld(i+nres)-vbldsc0(1,itype(i))
503             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
504      &       i,(dC(j,i+nres),j=1,3),xx
505           endif
506         enddo
507       endif
508 c Second correction (rotational lengthening)
509 c      do iter=1,maxiter
510       diffmax=0.0d0
511       ind=0
512       do i=nnt,nct-1
513         ind=ind+1
514         blen2 = scalar(dc(1,i),dc(1,i))
515         ppvec(ind)=2*vbl**2-blen2
516         diffbond=dabs(vbl-dsqrt(blen2))
517         if (diffbond.gt.diffmax) diffmax=diffbond
518         if (ppvec(ind).gt.0.0d0) then
519           ppvec(ind)=dsqrt(ppvec(ind))
520         else
521           ppvec(ind)=0.0d0
522         endif
523         if (lprn) then
524           write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
525         endif
526       enddo
527       do i=nnt,nct
528         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
529           ind=ind+1
530           blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
531           ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2
532           diffbond=dabs(vbldsc0(1,itype(i))-dsqrt(blen2))
533           if (diffbond.gt.diffmax) diffmax=diffbond
534           if (ppvec(ind).gt.0.0d0) then
535             ppvec(ind)=dsqrt(ppvec(ind))
536           else
537             ppvec(ind)=0.0d0
538           endif
539           if (lprn) then
540             write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
541           endif
542         endif
543       enddo
544       if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
545       if (diffmax.lt.difftol) goto 10
546       do i=1,dimen
547         Td(i)=0.0d0
548         do j=1,nbond
549           Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
550         enddo
551       enddo 
552       do i=1,dimen
553         zapas(i)=Td(i)
554         do j=1,dimen
555           zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
556         enddo
557       enddo
558       do j=1,3
559         dc(j,0)=zapas(j)
560         dc_work(j)=zapas(j)
561       enddo
562       ind=3
563       do i=nnt,nct-1
564         do j=1,3
565           dc(j,i)=zapas(ind+j)
566           dc_work(ind+j)=zapas(ind+j)
567         enddo
568         ind=ind+3
569       enddo
570       do i=nnt,nct
571         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
572           do j=1,3
573             dc(j,i+nres)=zapas(ind+j)
574             dc_work(ind+j)=zapas(ind+j)
575           enddo
576           ind=ind+3
577         endif
578       enddo 
579 c   Building the chain from the newly calculated coordinates    
580       call chainbuild_cart
581       if(ntwe.ne.0) then
582       if (large.and. mod(itime,ntwe).eq.0) then
583         write (iout,*) "Cartesian and internal coordinates: step 1"
584         call cartprint
585         call intout
586         write (iout,'(a)') "Potential forces"
587         do i=0,nres
588           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),
589      &    (-gxcart(j,i),j=1,3)
590         enddo
591         write (iout,'(a)') "Stochastic forces"
592         do i=0,nres
593           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),
594      &    (stochforc(j,i+nres),j=1,3)
595         enddo
596         write (iout,'(a)') "Velocities"
597         do i=0,nres
598           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
599      &    (d_t(j,i+nres),j=1,3)
600         enddo
601       endif
602       endif
603       if (lprn) then
604         write (iout,*) "After correction for rotational lengthening"
605         write (iout,*) "New coordinates",
606      &  " and differences between actual and standard bond lengths"
607         ind=0
608         do i=nnt,nct-1
609           ind=ind+1
610           xx=vbld(i+1)-vbl
611           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
612      &        i,(dC(j,i),j=1,3),xx
613         enddo
614         do i=nnt,nct
615           if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
616             ind=ind+1
617             xx=vbld(i+nres)-vbldsc0(1,itype(i))
618             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
619      &       i,(dC(j,i+nres),j=1,3),xx
620           endif
621         enddo
622       endif
623 c      ENDDO
624 c      write (iout,*) "Too many attempts at correcting the bonds"
625 c      stop
626    10 continue
627 #ifdef MPI
628       tt0 =MPI_Wtime()
629 #else
630       tt0 = tcpu()
631 #endif
632 c Calculate energy and forces
633       call zerograd
634       call etotal(potEcomp)
635       potE=potEcomp(0)-potEcomp(20)
636       call cartgrad
637       totT=totT+d_time
638 c  Calculate the kinetic and total energy and the kinetic temperature
639       call kinetic(EK)
640 #ifdef MPI
641       t_enegrad=t_enegrad+MPI_Wtime()-tt0
642 #else
643       t_enegrad=t_enegrad+tcpu()-tt0
644 #endif
645       totE=EK+potE
646       kinetic_T=2.0d0/(dimen*Rb)*EK
647       return
648       end
649 c-------------------------------------------------------------------------------
650       subroutine velverlet_step(itime)
651 c-------------------------------------------------------------------------------
652 c  Perform a single velocity Verlet step; the time step can be rescaled if 
653 c  increments in accelerations exceed the threshold
654 c-------------------------------------------------------------------------------
655       implicit real*8 (a-h,o-z)
656       include 'DIMENSIONS'
657 #ifdef MPI
658       include 'mpif.h'
659       integer ierror,ierrcode
660 #endif
661       include 'COMMON.CONTROL'
662       include 'COMMON.VAR'
663       include 'COMMON.MD'
664 #ifndef LANG0
665       include 'COMMON.LANGEVIN'
666 #else
667       include 'COMMON.LANGEVIN.lang0'
668 #endif
669       include 'COMMON.CHAIN'
670       include 'COMMON.DERIV'
671       include 'COMMON.GEO'
672       include 'COMMON.LOCAL'
673       include 'COMMON.INTERACT'
674       include 'COMMON.IOUNITS'
675       include 'COMMON.NAMES'
676       include 'COMMON.TIME1'
677       include 'COMMON.MUCA'
678       double precision vcm(3),incr(3)
679       double precision cm(3),L(3)
680       integer ilen,count,rstcount
681       external ilen
682       character*50 tytul
683       integer maxcount_scale /20/
684       common /gucio/ cm
685       double precision stochforcvec(MAXRES6)
686       common /stochcalc/ stochforcvec
687       integer itime
688       logical scale
689 c
690       scale=.true.
691       icount_scale=0
692       if (lang.eq.1) then
693         call sddir_precalc
694       else if (lang.eq.2 .or. lang.eq.3) then
695         call stochastic_force(stochforcvec)
696       endif
697       itime_scal=0
698       do while (scale)
699         icount_scale=icount_scale+1
700         if (icount_scale.gt.maxcount_scale) then
701           write (iout,*) 
702      &      "ERROR: too many attempts at scaling down the time step. ",
703      &      "amax=",amax,"epdrift=",epdrift,
704      &      "damax=",damax,"edriftmax=",edriftmax,
705      &      "d_time=",d_time
706           call flush(iout)
707 #ifdef MPI
708           call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
709 #endif
710           stop
711         endif
712 c First step of the velocity Verlet algorithm
713         if (lang.eq.2) then
714           call sd_verlet1
715         else if (lang.eq.3) then
716           call sd_verlet1_ciccotti
717         else if (lang.eq.1) then
718           call sddir_verlet1
719         else
720           call verlet1
721         endif     
722 c Build the chain from the newly calculated coordinates 
723         call chainbuild_cart
724         if (rattle) call rattle1
725         if (ntwe.ne.0) then
726         if (large.and. mod(itime,ntwe).eq.0) then
727           write (iout,*) "Cartesian and internal coordinates: step 1"
728           call cartprint
729           call intout
730           write (iout,*) "dC"
731           do i=0,nres
732             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
733      &      (dc(j,i+nres),j=1,3)
734           enddo
735           write (iout,*) "Accelerations"
736           do i=0,nres
737             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
738      &      (d_a(j,i+nres),j=1,3)
739           enddo
740           write (iout,*) "Velocities, step 1"
741           do i=0,nres
742             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
743      &      (d_t(j,i+nres),j=1,3)
744           enddo
745         endif
746         endif
747 #ifdef MPI
748         tt0 = MPI_Wtime()
749 #else
750         tt0 = tcpu()
751 #endif
752 c Calculate energy and forces
753         call zerograd
754         call etotal(potEcomp)
755         potE=potEcomp(0)-potEcomp(20)
756         call cartgrad
757 c Get the new accelerations
758         call lagrangian
759 #ifdef MPI
760         t_enegrad=t_enegrad+MPI_Wtime()-tt0
761 #else
762         t_enegrad=t_enegrad+tcpu()-tt0
763 #endif
764 c Determine maximum acceleration and scale down the timestep if needed
765         call max_accel
766         call predict_edrift(epdrift)
767         if (amax.gt.damax .or. epdrift.gt.edriftmax) then
768 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
769           scale=.true.
770           ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
771      &      /dlog(2.0d0)+1
772           itime_scal=itime_scal+ifac_time
773 c          fac_time=dmin1(damax/amax,0.5d0)
774           fac_time=0.5d0**ifac_time
775           d_time=d_time*fac_time
776           if (lang.eq.2 .or. lang.eq.3) then 
777 c            write (iout,*) "Calling sd_verlet_setup: 1"
778 c Rescale the stochastic forces and recalculate or restore 
779 c the matrices of tinker integrator
780             if (itime_scal.gt.maxflag_stoch) then
781               if (large) write (iout,'(a,i5,a)') 
782      &         "Calculate matrices for stochastic step;",
783      &         " itime_scal ",itime_scal
784               if (lang.eq.2) then
785                 call sd_verlet_p_setup
786               else
787                 call sd_verlet_ciccotti_setup
788               endif
789               write (iout,'(2a,i3,a,i3,1h.)') 
790      &         "Warning: cannot store matrices for stochastic",
791      &         " integration because the index",itime_scal,
792      &         " is greater than",maxflag_stoch
793               write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
794      &         " integration Langevin algorithm for better efficiency."
795             else if (flag_stoch(itime_scal)) then
796               if (large) write (iout,'(a,i5,a,l1)') 
797      &         "Restore matrices for stochastic step; itime_scal ",
798      &         itime_scal," flag ",flag_stoch(itime_scal)
799 #ifndef   LANG0
800               do i=1,dimen
801                 do j=1,dimen
802                   pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
803                   afric_mat(i,j)=afric0_mat(i,j,itime_scal)
804                   vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
805                   prand_mat(i,j)=prand0_mat(i,j,itime_scal)
806                   vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
807                   vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
808                 enddo
809               enddo
810 #endif
811             else
812               if (large) write (iout,'(2a,i5,a,l1)') 
813      &         "Calculate & store matrices for stochastic step;",
814      &         " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
815               if (lang.eq.2) then
816                 call sd_verlet_p_setup  
817               else
818                 call sd_verlet_ciccotti_setup
819               endif
820               flag_stoch(ifac_time)=.true.
821 #ifndef   LANG0
822               do i=1,dimen
823                 do j=1,dimen
824                   pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
825                   afric0_mat(i,j,itime_scal)=afric_mat(i,j)
826                   vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
827                   prand0_mat(i,j,itime_scal)=prand_mat(i,j)
828                   vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
829                   vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
830                 enddo
831               enddo
832 #endif
833             endif
834             fac_time=1.0d0/dsqrt(fac_time)
835             do i=1,dimen
836               stochforcvec(i)=fac_time*stochforcvec(i)
837             enddo
838           else if (lang.eq.1) then
839 c Rescale the accelerations due to stochastic forces
840             fac_time=1.0d0/dsqrt(fac_time)
841             do i=1,dimen
842               d_as_work(i)=d_as_work(i)*fac_time
843             enddo
844           endif
845           if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')  
846      &      "itime",itime," Timestep scaled down to ",
847      &      d_time," ifac_time",ifac_time," itime_scal",itime_scal
848         else 
849 c Second step of the velocity Verlet algorithm
850           if (lang.eq.2) then   
851             call sd_verlet2
852           else if (lang.eq.3) then
853             call sd_verlet2_ciccotti
854           else if (lang.eq.1) then
855             call sddir_verlet2
856           else
857             call verlet2
858           endif                     
859           if (rattle) call rattle2
860           totT=totT+d_time
861           if (d_time.ne.d_time0) then
862             d_time=d_time0
863             if (lang.eq.2 .or. lang.eq.3) then
864               if (large) write (iout,'(a)') 
865      &         "Restore original matrices for stochastic step"
866 c              write (iout,*) "Calling sd_verlet_setup: 2"
867 c Restore the matrices of tinker integrator if the time step has been restored
868 #ifndef   LANG0
869               do i=1,dimen
870                 do j=1,dimen
871                   pfric_mat(i,j)=pfric0_mat(i,j,0)
872                   afric_mat(i,j)=afric0_mat(i,j,0)
873                   vfric_mat(i,j)=vfric0_mat(i,j,0)
874                   prand_mat(i,j)=prand0_mat(i,j,0)
875                   vrand_mat1(i,j)=vrand0_mat1(i,j,0)
876                   vrand_mat2(i,j)=vrand0_mat2(i,j,0)
877                 enddo
878               enddo
879 #endif
880             endif
881           endif
882           scale=.false.
883         endif
884       enddo
885 c Calculate the kinetic and the total energy and the kinetic temperature
886       call kinetic(EK)
887       totE=EK+potE
888 c diagnostics
889 c      call kinetic1(EK1)
890 c      write (iout,*) "step",itime," EK",EK," EK1",EK1
891 c end diagnostics
892 c Couple the system to Berendsen bath if needed
893       if (tbf .and. lang.eq.0) then
894         call verlet_bath
895       endif
896       kinetic_T=2.0d0/(dimen*Rb)*EK
897 c Backup the coordinates, velocities, and accelerations
898       do i=0,2*nres
899         do j=1,3
900           dc_old(j,i)=dc(j,i)
901           d_t_old(j,i)=d_t(j,i)
902           d_a_old(j,i)=d_a(j,i)
903         enddo
904       enddo 
905       if (ntwe.ne.0) then
906       if (mod(itime,ntwe).eq.0 .and. large) then
907         write (iout,*) "Velocities, step 2"
908         do i=0,nres
909           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
910      &      (d_t(j,i+nres),j=1,3)
911         enddo
912       endif
913       endif
914       return
915       end
916 c-------------------------------------------------------------------------------
917       subroutine RESPA_step(itime)
918 c-------------------------------------------------------------------------------
919 c  Perform a single RESPA step.
920 c-------------------------------------------------------------------------------
921       implicit real*8 (a-h,o-z)
922       include 'DIMENSIONS'
923 #ifdef MPI
924       include 'mpif.h'
925 #endif
926       include 'COMMON.CONTROL'
927       include 'COMMON.VAR'
928       include 'COMMON.MD'
929 #ifndef LANG0
930       include 'COMMON.LANGEVIN'
931 #else
932       include 'COMMON.LANGEVIN.lang0'
933 #endif
934       include 'COMMON.CHAIN'
935       include 'COMMON.DERIV'
936       include 'COMMON.GEO'
937       include 'COMMON.LOCAL'
938       include 'COMMON.INTERACT'
939       include 'COMMON.IOUNITS'
940       include 'COMMON.NAMES'
941       include 'COMMON.TIME1'
942       double precision energia_short(0:n_ene),
943      & energia_long(0:n_ene)
944       double precision cm(3),L(3),vcm(3),incr(3)
945       integer ilen,count,rstcount
946       external ilen
947       character*50 tytul
948       integer maxcount_scale /10/
949       common /gucio/ cm
950       double precision stochforcvec(MAXRES6)
951       common /stochcalc/ stochforcvec
952       integer itime
953       logical scale
954       common /cipiszcze/ itt
955       itt=itime
956       large=(itime.gt.14600 .and. itime.lt.14700)
957       if (ntwe.ne.0) then
958       if (large.and. mod(itime,ntwe).eq.0) then
959         write (iout,*) "***************** RESPA itime",itime
960         write (iout,*) "Cartesian and internal coordinates: step 0"
961 c        call cartprint
962         call pdbout(0.0d0,"cipiszcze",iout)
963         call intout
964         write (iout,*) "Accelerations"
965         do i=0,nres
966           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
967      &      (d_a(j,i+nres),j=1,3)
968         enddo
969         write (iout,*) "Velocities, step 0"
970         do i=0,nres
971           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
972      &      (d_t(j,i+nres),j=1,3)
973         enddo
974       endif
975       endif
976 c
977 c Perform the initial RESPA step (increment velocities)
978 c      write (iout,*) "*********************** RESPA ini"
979       call RESPA_vel
980       if (ntwe.ne.0) then
981       if (mod(itime,ntwe).eq.0 .and. large) then
982         write (iout,*) "Velocities, end"
983         do i=0,nres
984           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
985      &      (d_t(j,i+nres),j=1,3)
986         enddo
987       endif
988       endif
989 c Compute the short-range forces
990 #ifdef MPI
991       tt0 =MPI_Wtime()
992 #else
993       tt0 = tcpu()
994 #endif
995       call zerograd
996       call etotal_short(energia_short)
997       if (large) write (iout,*) "energia_short",energia_short(0)
998       call cartgrad
999       call lagrangian
1000 #ifdef MPI
1001         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1002 #else
1003         t_enegrad=t_enegrad+tcpu()-tt0
1004 #endif
1005       do i=0,2*nres
1006         do j=1,3
1007           dc_old(j,i)=dc(j,i)
1008           d_t_old(j,i)=d_t(j,i)
1009           d_a_old(j,i)=d_a(j,i)
1010         enddo
1011       enddo 
1012       d_time0=d_time
1013 c Split the time step
1014       d_time=d_time/ntime_split 
1015 c Perform the short-range RESPSA steps (velocity Verlet increments of
1016 c positions and velocities using short-range forces)
1017 c      write (iout,*) "*********************** RESPA split"
1018       do itsplit=1,ntime_split
1019         if (lang.eq.1) then
1020           call sddir_precalc
1021         else if (lang.eq.2 .or. lang.eq.3) then
1022           call stochastic_force(stochforcvec)
1023         endif
1024 c First step of the velocity Verlet algorithm
1025         if (lang.eq.2) then
1026           call sd_verlet1
1027         else if (lang.eq.3) then
1028           call sd_verlet1_ciccotti
1029         else if (lang.eq.1) then
1030           call sddir_verlet1
1031         else
1032           call verlet1
1033         endif     
1034 c Build the chain from the newly calculated coordinates 
1035         call chainbuild_cart
1036         if (rattle) call rattle1
1037         if (ntwe.ne.0) then
1038         if (large.and. mod(itime,ntwe).eq.0) then
1039           write (iout,*) "Cartesian and internal coordinates: step 1"
1040           call pdbout(0.0d0,"cipiszcze",iout)
1041 c          call cartprint
1042           call intout
1043           write (iout,*) "Accelerations"
1044           do i=0,nres
1045             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1046      &        (d_a(j,i+nres),j=1,3)
1047           enddo
1048           write (iout,*) "Velocities, step 1"
1049           do i=0,nres
1050             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1051      &        (d_t(j,i+nres),j=1,3)
1052           enddo
1053         endif
1054         endif
1055 #ifdef MPI
1056         tt0 = MPI_Wtime()
1057 #else
1058         tt0 = tcpu()
1059 #endif
1060 c Calculate energy and forces
1061         call zerograd
1062         call etotal_short(energia_short)
1063       if (large) write (iout,*) "energia_short",energia_short(0)
1064         call cartgrad
1065 c Get the new accelerations
1066         call lagrangian
1067 #ifdef MPI
1068         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1069 #else
1070         t_enegrad=t_enegrad+tcpu()-tt0
1071 #endif
1072 c Second step of the velocity Verlet algorithm
1073         if (lang.eq.2) then     
1074           call sd_verlet2
1075         else if (lang.eq.3) then
1076           call sd_verlet2_ciccotti
1077         else if (lang.eq.1) then
1078           call sddir_verlet2
1079         else
1080           call verlet2
1081         endif               
1082         if (rattle) call rattle2
1083 c Backup the coordinates, velocities, and accelerations
1084         do i=0,2*nres
1085           do j=1,3
1086             dc_old(j,i)=dc(j,i)
1087             d_t_old(j,i)=d_t(j,i)
1088             d_a_old(j,i)=d_a(j,i)
1089           enddo
1090         enddo 
1091       enddo
1092 c Restore the time step
1093       d_time=d_time0
1094 c Compute long-range forces
1095 #ifdef MPI
1096       tt0 =MPI_Wtime()
1097 #else
1098       tt0 = tcpu()
1099 #endif
1100       call zerograd
1101       call etotal_long(energia_long)
1102       if (large) write (iout,*) "energia_long",energia_long(0)
1103       call cartgrad
1104       call lagrangian
1105 #ifdef MPI
1106         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1107 #else
1108         t_enegrad=t_enegrad+tcpu()-tt0
1109 #endif
1110 c Compute accelerations from long-range forces
1111       if (ntwe.ne.0) then
1112       if (large.and. mod(itime,ntwe).eq.0) then
1113         write (iout,*) "Cartesian and internal coordinates: step 2"
1114 c        call cartprint
1115         call pdbout(0.0d0,"cipiszcze",iout)
1116         call intout
1117         write (iout,*) "Accelerations"
1118         do i=0,nres
1119           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1120      &      (d_a(j,i+nres),j=1,3)
1121         enddo
1122         write (iout,*) "Velocities, step 2"
1123         do i=0,nres
1124           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1125      &      (d_t(j,i+nres),j=1,3)
1126         enddo
1127       endif
1128       endif
1129 c Compute the final RESPA step (increment velocities)
1130 c      write (iout,*) "*********************** RESPA fin"
1131       call RESPA_vel
1132 c Compute the complete potential energy
1133       do i=0,n_ene
1134         potEcomp(i)=energia_short(i)+energia_long(i)
1135       enddo
1136       potE=potEcomp(0)-potEcomp(20)
1137 c      potE=energia_short(0)+energia_long(0)
1138       totT=totT+d_time
1139 c Calculate the kinetic and the total energy and the kinetic temperature
1140       call kinetic(EK)
1141       totE=EK+potE
1142 c Couple the system to Berendsen bath if needed
1143       if (tbf .and. lang.eq.0) then
1144         call verlet_bath
1145       endif
1146       kinetic_T=2.0d0/(dimen*Rb)*EK
1147 c Backup the coordinates, velocities, and accelerations
1148       if (ntwe.ne.0) then
1149       if (mod(itime,ntwe).eq.0 .and. large) then
1150         write (iout,*) "Velocities, end"
1151         do i=0,nres
1152           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1153      &      (d_t(j,i+nres),j=1,3)
1154         enddo
1155       endif
1156       endif
1157       return
1158       end
1159 c---------------------------------------------------------------------
1160       subroutine RESPA_vel
1161 c  First and last RESPA step (incrementing velocities using long-range
1162 c  forces).
1163       implicit real*8 (a-h,o-z)
1164       include 'DIMENSIONS'
1165       include 'COMMON.CONTROL'
1166       include 'COMMON.VAR'
1167       include 'COMMON.MD'
1168       include 'COMMON.CHAIN'
1169       include 'COMMON.DERIV'
1170       include 'COMMON.GEO'
1171       include 'COMMON.LOCAL'
1172       include 'COMMON.INTERACT'
1173       include 'COMMON.IOUNITS'
1174       include 'COMMON.NAMES'
1175       do j=1,3
1176         d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1177       enddo
1178       do i=nnt,nct-1
1179         do j=1,3
1180           d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1181         enddo
1182       enddo
1183       do i=nnt,nct
1184         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1185           inres=i+nres
1186           do j=1,3
1187             d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1188           enddo
1189         endif
1190       enddo 
1191       return
1192       end
1193 c-----------------------------------------------------------------
1194       subroutine verlet1
1195 c Applying velocity Verlet algorithm - step 1 to coordinates
1196       implicit real*8 (a-h,o-z)
1197       include 'DIMENSIONS'
1198       include 'COMMON.CONTROL'
1199       include 'COMMON.VAR'
1200       include 'COMMON.MD'
1201       include 'COMMON.CHAIN'
1202       include 'COMMON.DERIV'
1203       include 'COMMON.GEO'
1204       include 'COMMON.LOCAL'
1205       include 'COMMON.INTERACT'
1206       include 'COMMON.IOUNITS'
1207       include 'COMMON.NAMES'
1208       double precision adt,adt2
1209         
1210       do j=1,3
1211         adt=d_a_old(j,0)*d_time
1212         adt2=0.5d0*adt
1213         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1214         d_t_new(j,0)=d_t_old(j,0)+adt2
1215         d_t(j,0)=d_t_old(j,0)+adt
1216       enddo
1217       do i=nnt,nct-1    
1218         do j=1,3    
1219           adt=d_a_old(j,i)*d_time
1220           adt2=0.5d0*adt
1221           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1222           d_t_new(j,i)=d_t_old(j,i)+adt2
1223           d_t(j,i)=d_t_old(j,i)+adt
1224         enddo
1225       enddo
1226       do i=nnt,nct
1227         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1228           inres=i+nres
1229           do j=1,3    
1230             adt=d_a_old(j,inres)*d_time
1231             adt2=0.5d0*adt
1232             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1233             d_t_new(j,inres)=d_t_old(j,inres)+adt2
1234             d_t(j,inres)=d_t_old(j,inres)+adt
1235           enddo
1236         endif      
1237       enddo 
1238       return
1239       end
1240 c---------------------------------------------------------------------
1241       subroutine verlet2
1242 c  Step 2 of the velocity Verlet algorithm: update velocities
1243       implicit real*8 (a-h,o-z)
1244       include 'DIMENSIONS'
1245       include 'COMMON.CONTROL'
1246       include 'COMMON.VAR'
1247       include 'COMMON.MD'
1248       include 'COMMON.CHAIN'
1249       include 'COMMON.DERIV'
1250       include 'COMMON.GEO'
1251       include 'COMMON.LOCAL'
1252       include 'COMMON.INTERACT'
1253       include 'COMMON.IOUNITS'
1254       include 'COMMON.NAMES'
1255       do j=1,3
1256         d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1257       enddo
1258       do i=nnt,nct-1
1259         do j=1,3
1260           d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1261         enddo
1262       enddo
1263       do i=nnt,nct
1264         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1265           inres=i+nres
1266           do j=1,3
1267             d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1268           enddo
1269         endif
1270       enddo 
1271       return
1272       end
1273 c-----------------------------------------------------------------
1274       subroutine sddir_precalc
1275 c Applying velocity Verlet algorithm - step 1 to coordinates        
1276       implicit real*8 (a-h,o-z)
1277       include 'DIMENSIONS'
1278       include 'COMMON.CONTROL'
1279       include 'COMMON.VAR'
1280       include 'COMMON.MD'
1281 #ifndef LANG0
1282       include 'COMMON.LANGEVIN'
1283 #else
1284       include 'COMMON.LANGEVIN.lang0'
1285 #endif
1286       include 'COMMON.CHAIN'
1287       include 'COMMON.DERIV'
1288       include 'COMMON.GEO'
1289       include 'COMMON.LOCAL'
1290       include 'COMMON.INTERACT'
1291       include 'COMMON.IOUNITS'
1292       include 'COMMON.NAMES'
1293       double precision stochforcvec(MAXRES6)
1294       common /stochcalc/ stochforcvec
1295 c
1296 c Compute friction and stochastic forces
1297 c
1298       call friction_force
1299       call stochastic_force(stochforcvec) 
1300 c
1301 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1302 c forces (d_as_work)
1303 c
1304 #ifdef OLD_GINV
1305       do i=1,dimen
1306         d_af_work(i)=0.0d0
1307         d_as_work(i)=0.0d0
1308         do j=1,dimen
1309           d_af_work(i)=d_af_work(i)+Ginv(i,j)*fric_work(j)
1310           d_as_work(i)=d_as_work(i)+Ginv(i,j)*stochforcvec(j)
1311         enddo
1312       enddo
1313 #else
1314         call ginv_mult(fric_work, d_af_work)
1315         call ginv_mult(stochforcvec, d_as_work)
1316 #endif
1317       return
1318       end
1319 c---------------------------------------------------------------------
1320       subroutine sddir_verlet1
1321 c Applying velocity Verlet algorithm - step 1 to velocities        
1322       implicit real*8 (a-h,o-z)
1323       include 'DIMENSIONS'
1324       include 'COMMON.CONTROL'
1325       include 'COMMON.VAR'
1326       include 'COMMON.MD'
1327 #ifndef LANG0
1328       include 'COMMON.LANGEVIN'
1329 #else
1330       include 'COMMON.LANGEVIN.lang0'
1331 #endif
1332       include 'COMMON.CHAIN'
1333       include 'COMMON.DERIV'
1334       include 'COMMON.GEO'
1335       include 'COMMON.LOCAL'
1336       include 'COMMON.INTERACT'
1337       include 'COMMON.IOUNITS'
1338       include 'COMMON.NAMES'
1339 c Revised 3/31/05 AL: correlation between random contributions to 
1340 c position and velocity increments included.
1341       double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1342       double precision adt,adt2
1343 c
1344 c Add the contribution from BOTH friction and stochastic force to the
1345 c coordinates, but ONLY the contribution from the friction forces to velocities
1346 c
1347       do j=1,3
1348         adt=(d_a_old(j,0)+d_af_work(j))*d_time
1349         adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1350         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1351         d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1352         d_t(j,0)=d_t_old(j,0)+adt
1353       enddo
1354       ind=3
1355       do i=nnt,nct-1    
1356         do j=1,3    
1357           adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1358           adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1359           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1360           d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1361           d_t(j,i)=d_t_old(j,i)+adt
1362         enddo
1363         ind=ind+3
1364       enddo
1365       do i=nnt,nct
1366         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1367           inres=i+nres
1368           do j=1,3    
1369             adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1370             adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1371             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1372             d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1373             d_t(j,inres)=d_t_old(j,inres)+adt
1374           enddo
1375           ind=ind+3
1376         endif      
1377       enddo 
1378       return
1379       end
1380 c---------------------------------------------------------------------
1381       subroutine sddir_verlet2
1382 c  Calculating the adjusted velocities for accelerations
1383       implicit real*8 (a-h,o-z)
1384       include 'DIMENSIONS'
1385       include 'COMMON.CONTROL'
1386       include 'COMMON.VAR'
1387       include 'COMMON.MD'
1388 #ifndef LANG0
1389       include 'COMMON.LANGEVIN'
1390 #else
1391       include 'COMMON.LANGEVIN.lang0'
1392 #endif
1393       include 'COMMON.CHAIN'
1394       include 'COMMON.DERIV'
1395       include 'COMMON.GEO'
1396       include 'COMMON.LOCAL'
1397       include 'COMMON.INTERACT'
1398       include 'COMMON.IOUNITS'
1399       include 'COMMON.NAMES'
1400       double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1401       double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1402 c Revised 3/31/05 AL: correlation between random contributions to 
1403 c position and velocity increments included.
1404 c The correlation coefficients are calculated at low-friction limit.
1405 c Also, friction forces are now not calculated with new velocities.
1406
1407 c      call friction_force
1408       call stochastic_force(stochforcvec) 
1409 c
1410 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1411 c forces (d_as_work)
1412 c
1413 #ifdef OLD_GINV
1414       do i=1,dimen
1415 c        d_af_work(i)=0.0d0
1416         d_as_work1(i)=0.0d0
1417         do j=1,dimen
1418 c          d_af_work(i)=d_af_work(i)+Ginv(i,j)*fric_work(j)
1419           d_as_work1(i)=d_as_work1(i)+Ginv(i,j)*stochforcvec(j)
1420         enddo
1421       enddo
1422 #else
1423         call ginv_mult(stochforcvec, d_as_work1)
1424 #endif
1425
1426 c
1427 c Update velocities
1428 c
1429       do j=1,3
1430         d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1431      &    +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1432       enddo
1433       ind=3
1434       do i=nnt,nct-1
1435         do j=1,3
1436           d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1437      &     +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1438         enddo
1439         ind=ind+3
1440       enddo
1441       do i=nnt,nct
1442         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1443           inres=i+nres
1444           do j=1,3
1445             d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1446      &       +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1447      &       +cos60*d_as_work1(ind+j))*d_time
1448           enddo
1449           ind=ind+3
1450         endif
1451       enddo 
1452       return
1453       end
1454 c---------------------------------------------------------------------
1455       subroutine max_accel
1456 c
1457 c Find the maximum difference in the accelerations of the the sites
1458 c at the beginning and the end of the time step.
1459 c
1460       implicit real*8 (a-h,o-z)
1461       include 'DIMENSIONS'
1462       include 'COMMON.CONTROL'
1463       include 'COMMON.VAR'
1464       include 'COMMON.MD'
1465       include 'COMMON.CHAIN'
1466       include 'COMMON.DERIV'
1467       include 'COMMON.GEO'
1468       include 'COMMON.LOCAL'
1469       include 'COMMON.INTERACT'
1470       include 'COMMON.IOUNITS'
1471       double precision aux(3),accel(3)
1472       do j=1,3
1473         aux(j)=d_a(j,0)-d_a_old(j,0)
1474       enddo 
1475       amax=0.0d0
1476       do i=nnt,nct
1477 c Backbone
1478         if (i.lt.nct) then
1479           do j=1,3
1480             accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1481             if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1482           enddo
1483         endif
1484 c Side chains
1485         do j=1,3
1486           accel(j)=aux(j)
1487         enddo
1488         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1489           do j=1,3 
1490             accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1491           enddo
1492         endif
1493         do j=1,3
1494           if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1495         enddo
1496         do j=1,3
1497           aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1498         enddo
1499       enddo
1500       return
1501       end       
1502 c---------------------------------------------------------------------
1503       subroutine predict_edrift(epdrift)
1504 c
1505 c Predict the drift of the potential energy
1506 c
1507       implicit real*8 (a-h,o-z)
1508       include 'DIMENSIONS'
1509       include 'COMMON.CONTROL'
1510       include 'COMMON.VAR'
1511       include 'COMMON.MD'
1512       include 'COMMON.CHAIN'
1513       include 'COMMON.DERIV'
1514       include 'COMMON.GEO'
1515       include 'COMMON.LOCAL'
1516       include 'COMMON.INTERACT'
1517       include 'COMMON.IOUNITS'
1518       include 'COMMON.MUCA'
1519       double precision epdrift,epdriftij
1520 c Drift of the potential energy
1521       epdrift=0.0d0
1522       do i=nnt,nct
1523 c Backbone
1524         if (i.lt.nct) then
1525           do j=1,3
1526             epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1527             if (lmuca) epdriftij=epdriftij*factor
1528 c            write (iout,*) "back",i,j,epdriftij
1529             if (epdriftij.gt.epdrift) epdrift=epdriftij 
1530           enddo
1531         endif
1532 c Side chains
1533         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1534           do j=1,3 
1535             epdriftij=
1536      &       dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1537             if (lmuca) epdriftij=epdriftij*factor
1538 c            write (iout,*) "side",i,j,epdriftij
1539             if (epdriftij.gt.epdrift) epdrift=epdriftij
1540           enddo
1541         endif
1542       enddo
1543       epdrift=0.5d0*epdrift*d_time*d_time
1544 c      write (iout,*) "epdrift",epdrift
1545       return
1546       end       
1547 c-----------------------------------------------------------------------
1548       subroutine verlet_bath
1549 c
1550 c  Coupling to the thermostat by using the Berendsen algorithm
1551 c
1552       implicit real*8 (a-h,o-z)
1553       include 'DIMENSIONS'
1554       include 'COMMON.CONTROL'
1555       include 'COMMON.VAR'
1556       include 'COMMON.MD'
1557       include 'COMMON.CHAIN'
1558       include 'COMMON.DERIV'
1559       include 'COMMON.GEO'
1560       include 'COMMON.LOCAL'
1561       include 'COMMON.INTERACT'
1562       include 'COMMON.IOUNITS'
1563       include 'COMMON.NAMES'
1564       double precision T_half,fact
1565
1566       T_half=2.0d0/(dimen*Rb)*EK
1567       fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1568 c      write(iout,*) "T_half", T_half
1569 c      write(iout,*) "EK", EK
1570 c      write(iout,*) "fact", fact                               
1571       do j=1,3
1572         d_t(j,0)=fact*d_t(j,0)
1573       enddo
1574       do i=nnt,nct-1
1575         do j=1,3
1576           d_t(j,i)=fact*d_t(j,i)
1577         enddo
1578       enddo
1579       do i=nnt,nct
1580         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1581           inres=i+nres
1582           do j=1,3
1583             d_t(j,inres)=fact*d_t(j,inres)
1584           enddo
1585         endif
1586       enddo 
1587       return
1588       end
1589 c---------------------------------------------------------
1590       subroutine init_MD
1591 c  Set up the initial conditions of a MD simulation
1592       implicit real*8 (a-h,o-z)
1593       include 'DIMENSIONS'
1594 #ifdef MP
1595       include 'mpif.h'
1596       include 'COMMON.SETUP'
1597       character*16 form
1598 #endif
1599       include 'COMMON.CONTROL'
1600       include 'COMMON.VAR'
1601       include 'COMMON.MD'
1602 #ifndef LANG0
1603       include 'COMMON.LANGEVIN'
1604 #else
1605       include 'COMMON.LANGEVIN.lang0'
1606 #endif
1607       include 'COMMON.CHAIN'
1608       include 'COMMON.DERIV'
1609       include 'COMMON.GEO'
1610       include 'COMMON.LOCAL'
1611       include 'COMMON.INTERACT'
1612       include 'COMMON.IOUNITS'
1613       include 'COMMON.NAMES'
1614       include 'COMMON.REMD'
1615       real*8 energia_long(0:n_ene),
1616      &  energia_short(0:n_ene),vcm(3),incr(3)
1617       double precision cm(3),L(3),xv,sigv,lowb,highb
1618       double precision varia(maxvar)
1619       character*256 qstr
1620       integer ilen
1621       external ilen
1622       character*50 tytul
1623       logical file_exist
1624       common /gucio/ cm
1625       d_time0=d_time
1626 c      write(iout,*) "d_time", d_time
1627 c Compute the standard deviations of stochastic forces for Langevin dynamics
1628 c if the friction coefficients do not depend on surface area
1629       if (lang.gt.0 .and. .not.surfarea) then
1630         do i=nnt,nct-1
1631           stdforcp(i)=stdfp*dsqrt(gamp)
1632         enddo
1633         do i=nnt,nct
1634           stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(iabs(itype(i))))
1635         enddo
1636       endif
1637 c Open the pdb file for snapshotshots
1638 #ifdef MPI
1639       if(mdpdb) then
1640         if (ilen(tmpdir).gt.0) 
1641      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1642      &      liczba(:ilen(liczba))//".pdb")
1643         open(ipdb,
1644      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1645      &  //".pdb")
1646       else
1647 #ifdef NOXDR
1648         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) 
1649      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1650      &      liczba(:ilen(liczba))//".x")
1651         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1652      &  //".x"
1653 #else
1654         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) 
1655      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1656      &      liczba(:ilen(liczba))//".cx")
1657         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1658      &  //".cx"
1659 #endif
1660       endif
1661 #else
1662       if(mdpdb) then
1663          if (ilen(tmpdir).gt.0) 
1664      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1665          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1666       else
1667          if (ilen(tmpdir).gt.0) 
1668      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1669          cartname=prefix(:ilen(prefix))//"_MD.cx"
1670       endif
1671 #endif
1672       if (usampl) then
1673         write (qstr,'(256(1h ))')
1674         ipos=1
1675         do i=1,nfrag
1676           iq = qinfrag(i,iset)*10
1677           iw = wfrag(i,iset)/100
1678           if (iw.gt.0) then
1679             if(me.eq.king.or..not.out1file)
1680      &       write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1681             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1682             ipos=ipos+7
1683           endif
1684         enddo
1685         do i=1,npair
1686           iq = qinpair(i,iset)*10
1687           iw = wpair(i,iset)/100
1688           if (iw.gt.0) then
1689             if(me.eq.king.or..not.out1file)
1690      &       write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1691             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1692             ipos=ipos+7
1693           endif
1694         enddo
1695 c        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1696 #ifdef NOXDR
1697 c        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1698 #else
1699 c        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1700 #endif
1701 c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1702       endif
1703       icg=1
1704       if (rest) then
1705        if (restart1file) then
1706          if (me.eq.king)
1707      &     inquire(file=mremd_rst_name,exist=file_exist)
1708            write (*,*) me," Before broadcast: file_exist",file_exist
1709          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1710      &          IERR)
1711          write (*,*) me," After broadcast: file_exist",file_exist
1712 c        inquire(file=mremd_rst_name,exist=file_exist)
1713         if(me.eq.king.or..not.out1file)
1714      &   write(iout,*) "Initial state read by master and distributed"
1715        else
1716          if (ilen(tmpdir).gt.0)
1717      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1718      &      //liczba(:ilen(liczba))//'.rst')
1719         inquire(file=rest2name,exist=file_exist)
1720        endif
1721        if(file_exist) then
1722          if(.not.restart1file) then
1723            if(me.eq.king.or..not.out1file)
1724      &      write(iout,*) "Initial state will be read from file ",
1725      &      rest2name(:ilen(rest2name))
1726            call readrst
1727          endif  
1728          call rescale_weights(t_bath)
1729        else
1730         if(me.eq.king.or..not.out1file)then
1731          if (restart1file) then
1732           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1733      &       " does not exist"
1734          else
1735           write(iout,*) "File ",rest2name(:ilen(rest2name)),
1736      &       " does not exist"
1737          endif
1738          write(iout,*) "Initial velocities randomly generated"
1739         endif
1740         call random_vel
1741         totT=0.0d0
1742        endif
1743       else
1744 c Generate initial velocities
1745         if(me.eq.king.or..not.out1file)
1746      &   write(iout,*) "Initial velocities randomly generated"
1747         call random_vel
1748         totT=0.0d0
1749       endif
1750 c      rest2name = prefix(:ilen(prefix))//'.rst'
1751       if(me.eq.king.or..not.out1file)then
1752        write(iout,*) "Initial backbone velocities"
1753        do i=nnt,nct-1
1754         write(iout,*) (d_t(j,i),j=1,3)
1755        enddo
1756        write(iout,*) "Initial side-chain velocities"
1757        do i=nnt,nct
1758         write(iout,*) (d_t(j,i+nres),j=1,3)
1759        enddo                     
1760 c  Zeroing the total angular momentum of the system
1761        write(iout,*) "Calling the zero-angular 
1762      & momentum subroutine"
1763       endif
1764       call inertia_tensor  
1765 c  Getting the potential energy and forces and velocities and accelerations
1766       call vcm_vel(vcm)
1767 c      write (iout,*) "velocity of the center of the mass:"
1768 c      write (iout,*) (vcm(j),j=1,3)
1769       do j=1,3
1770         d_t(j,0)=d_t(j,0)-vcm(j)
1771       enddo
1772 c Removing the velocity of the center of mass
1773       call vcm_vel(vcm)
1774       if(me.eq.king.or..not.out1file)then
1775        write (iout,*) "vcm right after adjustment:"
1776        write (iout,*) (vcm(j),j=1,3) 
1777       endif
1778       if (.not.rest) then               
1779          call chainbuild
1780          if(iranconf.ne.0) then
1781           if (overlapsc) then 
1782            print *, 'Calling OVERLAP_SC'
1783            call overlap_sc(fail)
1784           endif 
1785
1786           if (searchsc) then 
1787            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1788            print *,'SC_move',nft_sc,etot
1789            if(me.eq.king.or..not.out1file)
1790      &      write(iout,*) 'SC_move',nft_sc,etot
1791           endif 
1792
1793           if(dccart)then
1794            print *, 'Calling MINIM_DC'
1795            call minim_dc(etot,iretcode,nfun)
1796           else
1797            call geom_to_var(nvar,varia)
1798            print *,'Calling MINIMIZE.'
1799            call minimize(etot,varia,iretcode,nfun)
1800            call var_to_geom(nvar,varia)
1801           endif
1802           if(me.eq.king.or..not.out1file)
1803      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1804          endif
1805       endif       
1806       call chainbuild_cart
1807       call kinetic(EK)
1808       if (tbf) then
1809         call verlet_bath(EK)
1810       endif       
1811       kinetic_T=2.0d0/(dimen*Rb)*EK
1812       if(me.eq.king.or..not.out1file)then
1813        call cartprint
1814        call intout
1815       endif
1816       call zerograd
1817       call etotal(potEcomp)
1818       potE=potEcomp(0)
1819       call cartgrad
1820       call lagrangian
1821       call max_accel
1822       if (amax*d_time .gt. dvmax) d_time=d_time*dvmax/amax
1823       if(me.eq.king.or..not.out1file)then 
1824        write(iout,*) "Potential energy and its components"
1825        write(iout,*) (potEcomp(i),i=0,n_ene)
1826       endif
1827       potE=potEcomp(0)-potEcomp(20)
1828       totE=EK+potE
1829       itime=0
1830       if (ntwe.ne.0) call statout(itime)
1831       if(me.eq.king.or..not.out1file)
1832      &  write (iout,*) "Initial:",
1833      &   " Kinetic energy",EK," potential energy",potE, 
1834      &   " total energy",totE," maximum acceleration ",
1835      &   amax
1836       if (large) then
1837         write (iout,*) "Initial coordinates"
1838         do i=1,nres
1839           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1840      &    (c(j,i+nres),j=1,3)
1841         enddo
1842         write (iout,*) "Initial dC"
1843         do i=0,nres
1844           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1845      &    (dc(j,i+nres),j=1,3)
1846         enddo
1847         write (iout,*) "Initial velocities"
1848         do i=0,nres
1849           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1850      &    (d_t(j,i+nres),j=1,3)
1851         enddo
1852         write (iout,*) "Initial accelerations"
1853         do i=0,nres
1854           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1855      &    (d_a(j,i+nres),j=1,3)
1856         enddo
1857       endif
1858       do i=0,2*nres
1859         do j=1,3
1860           dc_old(j,i)=dc(j,i)
1861           d_t_old(j,i)=d_t(j,i)
1862           d_a_old(j,i)=d_a(j,i)
1863         enddo
1864 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1865       enddo 
1866       if (RESPA) then
1867 #ifdef MPI
1868       tt0 =MPI_Wtime()
1869 #else
1870       tt0 = tcpu()
1871 #endif
1872         call zerograd
1873         call etotal_long(energia_long)
1874       if (large) write (iout,*) "energia_long",energia_long(0)
1875         call cartgrad
1876         call lagrangian
1877 #ifdef MPI
1878         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1879 #else
1880         t_enegrad=t_enegrad+tcpu()-tt0
1881 #endif
1882 c        call etotal_short(energia_short)
1883 c        write (iout,*) "energia_long",energia_long(0),
1884 c     &   " energia_short",energia_short(0),
1885 c     &   " total",energia_long(0)+energia_short(0)
1886       endif
1887       return
1888       end
1889 c-----------------------------------------------------------
1890       subroutine random_vel
1891       implicit real*8 (a-h,o-z)
1892       include 'DIMENSIONS'
1893       include 'COMMON.CONTROL'
1894       include 'COMMON.VAR'
1895       include 'COMMON.MD'
1896 #ifndef LANG0
1897       include 'COMMON.LANGEVIN'
1898 #else
1899       include 'COMMON.LANGEVIN.lang0'
1900 #endif
1901       include 'COMMON.CHAIN'
1902       include 'COMMON.DERIV'
1903       include 'COMMON.GEO'
1904       include 'COMMON.LOCAL'
1905       include 'COMMON.INTERACT'
1906       include 'COMMON.IOUNITS'
1907       include 'COMMON.NAMES'
1908       include 'COMMON.TIME1'
1909       double precision xv,sigv,lowb,highb
1910 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m 
1911 c First generate velocities in the eigenspace of the G matrix
1912 c      write (iout,*) "Calling random_vel"
1913       xv=0.0d0
1914       do i=1,dimen
1915         sigv=dsqrt((Rb*t_bath)/geigen(i))
1916         lowb=-5*sigv
1917         highb=5*sigv
1918         d_t_work_new(i)=anorm_distr(xv,sigv,lowb,highb)
1919       enddo
1920 c      Ek1=0.0d0
1921 c      do i=1,dimen
1922 c        Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(i)**2
1923 c      enddo
1924 c Transform velocities to UNRES coordinate space
1925       do i=1,dimen
1926         d_t_work(i)=0.0d0
1927         do j=1,dimen
1928           d_t_work(i)=d_t_work(i)+Gvec(i,j)*d_t_work_new(j)
1929         enddo
1930       enddo
1931 c Transfer to the d_t vector
1932       do j=1,3
1933         d_t(j,0)=d_t_work(j)
1934       enddo 
1935       ind=3
1936       do i=nnt,nct-1
1937         do j=1,3 
1938           ind=ind+1
1939           if (itype(i).ne.ntyp1 .and. itype(i+1).ne.ntyp1) then
1940             d_t(j,i)=d_t_work(ind)
1941           else
1942             d_t(j,i)=0.0d0
1943           endif
1944         enddo
1945       enddo
1946       do i=nnt,nct
1947         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1948           do j=1,3
1949             ind=ind+1
1950             d_t(j,i+nres)=d_t_work(ind)
1951           enddo
1952         endif
1953       enddo
1954 c      call kinetic(EK)
1955 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1956 c     &  2.0d0/(dimen*Rb)*EK,2.0d0/(dimen*Rb)*EK1
1957       return
1958       end
1959 c-----------------------------------------------------------
1960       subroutine sd_verlet_p_setup
1961 c Sets up the parameters of stochastic Verlet algorithm       
1962       implicit real*8 (a-h,o-z)
1963       include 'DIMENSIONS'
1964 #ifdef MPI
1965       include 'mpif.h'
1966 #endif
1967       include 'COMMON.CONTROL'
1968       include 'COMMON.VAR'
1969       include 'COMMON.MD'
1970 #ifndef LANG0
1971       include 'COMMON.LANGEVIN'
1972 #else
1973       include 'COMMON.LANGEVIN.lang0'
1974 #endif
1975       include 'COMMON.CHAIN'
1976       include 'COMMON.DERIV'
1977       include 'COMMON.GEO'
1978       include 'COMMON.LOCAL'
1979       include 'COMMON.INTERACT'
1980       include 'COMMON.IOUNITS'
1981       include 'COMMON.NAMES'
1982       include 'COMMON.TIME1'
1983       double precision emgdt(MAXRES6),
1984      & pterm,vterm,rho,rhoc,vsig,
1985      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1986      & afric_vec(MAXRES6),prand_vec(MAXRES6),
1987      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1988       logical lprn /.false./
1989       double precision zero /1.0d-8/, gdt_radius /0.05d0/ 
1990       double precision ktm
1991 #ifdef MPI
1992       tt0 = MPI_Wtime()
1993 #else
1994       tt0 = tcpu()
1995 #endif
1996 c
1997 c AL 8/17/04 Code adapted from tinker
1998 c
1999 c Get the frictional and random terms for stochastic dynamics in the
2000 c eigenspace of mass-scaled UNRES friction matrix
2001 c
2002       do i = 1, dimen
2003             gdt = fricgam(i) * d_time
2004 c
2005 c Stochastic dynamics reduces to simple MD for zero friction
2006 c
2007             if (gdt .le. zero) then
2008                pfric_vec(i) = 1.0d0
2009                vfric_vec(i) = d_time
2010                afric_vec(i) = 0.5d0 * d_time * d_time
2011                prand_vec(i) = 0.0d0
2012                vrand_vec1(i) = 0.0d0
2013                vrand_vec2(i) = 0.0d0
2014 c
2015 c Analytical expressions when friction coefficient is large
2016 c
2017             else 
2018                if (gdt .ge. gdt_radius) then
2019                   egdt = dexp(-gdt)
2020                   pfric_vec(i) = egdt
2021                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2022                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2023                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2024                   vterm = 1.0d0 - egdt**2
2025                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2026 c
2027 c Use series expansions when friction coefficient is small
2028 c
2029                else
2030                   gdt2 = gdt * gdt
2031                   gdt3 = gdt * gdt2
2032                   gdt4 = gdt2 * gdt2
2033                   gdt5 = gdt2 * gdt3
2034                   gdt6 = gdt3 * gdt3
2035                   gdt7 = gdt3 * gdt4
2036                   gdt8 = gdt4 * gdt4
2037                   gdt9 = gdt4 * gdt5
2038                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2039      &                          - gdt5/120.0d0 + gdt6/720.0d0
2040      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
2041      &                          - gdt9/362880.0d0) / fricgam(i)**2
2042                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2043                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2044                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2045      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2046      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2047      &                       + 127.0d0*gdt9/90720.0d0
2048                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2049      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2050      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2051      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2052                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2053      &                       - 17.0d0*gdt2/1280.0d0
2054      &                       + 17.0d0*gdt3/6144.0d0
2055      &                       + 40967.0d0*gdt4/34406400.0d0
2056      &                       - 57203.0d0*gdt5/275251200.0d0
2057      &                       - 1429487.0d0*gdt6/13212057600.0d0)
2058                end if
2059 c
2060 c Compute the scaling factors of random terms for the nonzero friction case
2061 c
2062                ktm = 0.5d0*d_time/fricgam(i)
2063                psig = dsqrt(ktm*pterm) / fricgam(i)
2064                vsig = dsqrt(ktm*vterm)
2065                rhoc = dsqrt(1.0d0 - rho*rho)
2066                prand_vec(i) = psig 
2067                vrand_vec1(i) = vsig * rho 
2068                vrand_vec2(i) = vsig * rhoc
2069             end if
2070       end do
2071       if (lprn) then
2072       write (iout,*) 
2073      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2074      &  " vrand_vec2"
2075       do i=1,dimen
2076         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2077      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2078       enddo
2079       endif
2080 c
2081 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2082 c
2083 #ifndef   LANG0
2084       call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
2085       call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
2086       call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
2087       call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
2088       call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec1,vrand_mat1)
2089       call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
2090 #endif
2091 c      call eigtransf1(dimen,maxres6,mt3mt2,pfric_vec,pfric_mat)
2092 c      call eigtransf1(dimen,maxres6,mt3mt2,vfric_vec,vfric_mat)
2093 c      call eigtransf1(dimen,maxres6,mt3mt2,afric_vec,afric_mat)
2094 c      call eigtransf1(dimen,maxres6,mt3mt1,prand_vec,prand_mat)
2095 c      call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec1,vrand_mat1)
2096 c      call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec2,vrand_mat2)
2097 #ifdef MPI
2098       t_sdsetup=t_sdsetup+MPI_Wtime()
2099 #else
2100       t_sdsetup=t_sdsetup+tcpu()-tt0
2101 #endif
2102       return
2103       end
2104 c-------------------------------------------------------------      
2105       subroutine eigtransf1(n,ndim,ab,d,c)
2106       implicit none
2107       integer n,ndim
2108       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2109       integer i,j,k
2110       do i=1,n
2111         do j=1,n
2112           c(i,j)=0.0d0
2113           do k=1,n
2114             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2115           enddo
2116         enddo
2117       enddo
2118       return
2119       end
2120 c-------------------------------------------------------------      
2121       subroutine eigtransf(n,ndim,a,b,d,c)
2122       implicit none
2123       integer n,ndim
2124       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2125       integer i,j,k
2126       do i=1,n
2127         do j=1,n
2128           c(i,j)=0.0d0
2129           do k=1,n
2130             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2131           enddo
2132         enddo
2133       enddo
2134       return
2135       end
2136 c-------------------------------------------------------------      
2137       subroutine sd_verlet1
2138 c Applying stochastic velocity Verlet algorithm - step 1 to velocities        
2139       implicit real*8 (a-h,o-z)
2140       include 'DIMENSIONS'
2141       include 'COMMON.CONTROL'
2142       include 'COMMON.VAR'
2143       include 'COMMON.MD'
2144 #ifndef LANG0
2145       include 'COMMON.LANGEVIN'
2146 #else
2147       include 'COMMON.LANGEVIN.lang0'
2148 #endif
2149       include 'COMMON.CHAIN'
2150       include 'COMMON.DERIV'
2151       include 'COMMON.GEO'
2152       include 'COMMON.LOCAL'
2153       include 'COMMON.INTERACT'
2154       include 'COMMON.IOUNITS'
2155       include 'COMMON.NAMES'
2156       double precision stochforcvec(MAXRES6)
2157       common /stochcalc/ stochforcvec
2158       logical lprn /.false./
2159
2160 c      write (iout,*) "dc_old"
2161 c      do i=0,nres
2162 c        write (iout,'(i5,3f10.5,5x,3f10.5)') 
2163 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2164 c      enddo
2165       do j=1,3
2166         dc_work(j)=dc_old(j,0)
2167         d_t_work(j)=d_t_old(j,0)
2168         d_a_work(j)=d_a_old(j,0)
2169       enddo
2170       ind=3
2171       do i=nnt,nct-1
2172         do j=1,3
2173           dc_work(ind+j)=dc_old(j,i)
2174           d_t_work(ind+j)=d_t_old(j,i)
2175           d_a_work(ind+j)=d_a_old(j,i)
2176         enddo
2177         ind=ind+3
2178       enddo
2179       do i=nnt,nct
2180         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2181           do j=1,3
2182             dc_work(ind+j)=dc_old(j,i+nres)
2183             d_t_work(ind+j)=d_t_old(j,i+nres)
2184             d_a_work(ind+j)=d_a_old(j,i+nres)
2185           enddo
2186           ind=ind+3
2187         endif
2188       enddo
2189 #ifndef LANG0
2190       if (lprn) then
2191       write (iout,*) 
2192      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2193      &  " vrand_mat2"
2194       do i=1,dimen
2195         do j=1,dimen
2196           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2197      &      vfric_mat(i,j),afric_mat(i,j),
2198      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2199         enddo
2200       enddo
2201       endif
2202       do i=1,dimen
2203         ddt1=0.0d0
2204         ddt2=0.0d0
2205         do j=1,dimen
2206           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2207      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2208           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2209           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2210         enddo
2211         d_t_work_new(i)=ddt1+0.5d0*ddt2
2212         d_t_work(i)=ddt1+ddt2
2213       enddo
2214 #endif
2215       do j=1,3
2216         dc(j,0)=dc_work(j)
2217         d_t(j,0)=d_t_work(j)
2218       enddo
2219       ind=3     
2220       do i=nnt,nct-1    
2221         do j=1,3
2222           dc(j,i)=dc_work(ind+j)
2223           d_t(j,i)=d_t_work(ind+j)
2224         enddo
2225         ind=ind+3
2226       enddo
2227       do i=nnt,nct
2228         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2229           inres=i+nres
2230           do j=1,3
2231             dc(j,inres)=dc_work(ind+j)
2232             d_t(j,inres)=d_t_work(ind+j)
2233           enddo
2234           ind=ind+3
2235         endif      
2236       enddo 
2237       return
2238       end
2239 c--------------------------------------------------------------------------
2240       subroutine sd_verlet2
2241 c  Calculating the adjusted velocities for accelerations
2242       implicit real*8 (a-h,o-z)
2243       include 'DIMENSIONS'
2244       include 'COMMON.CONTROL'
2245       include 'COMMON.VAR'
2246       include 'COMMON.MD'
2247 #ifndef LANG0
2248       include 'COMMON.LANGEVIN'
2249 #else
2250       include 'COMMON.LANGEVIN.lang0'
2251 #endif
2252       include 'COMMON.CHAIN'
2253       include 'COMMON.DERIV'
2254       include 'COMMON.GEO'
2255       include 'COMMON.LOCAL'
2256       include 'COMMON.INTERACT'
2257       include 'COMMON.IOUNITS'
2258       include 'COMMON.NAMES'
2259       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2260       common /stochcalc/ stochforcvec
2261 c
2262 c Compute the stochastic forces which contribute to velocity change
2263 c
2264       call stochastic_force(stochforcvecV)
2265
2266 #ifndef LANG0
2267       do i=1,dimen
2268         ddt1=0.0d0
2269         ddt2=0.0d0
2270         do j=1,dimen
2271           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2272           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2273      &     vrand_mat2(i,j)*stochforcvecV(j)
2274         enddo
2275         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2276       enddo
2277 #endif
2278       do j=1,3
2279         d_t(j,0)=d_t_work(j)
2280       enddo
2281       ind=3
2282       do i=nnt,nct-1
2283         do j=1,3
2284           d_t(j,i)=d_t_work(ind+j)
2285         enddo
2286         ind=ind+3
2287       enddo
2288       do i=nnt,nct
2289         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2290           inres=i+nres
2291           do j=1,3
2292             d_t(j,inres)=d_t_work(ind+j)
2293           enddo
2294           ind=ind+3
2295         endif
2296       enddo 
2297       return
2298       end
2299 c-----------------------------------------------------------
2300       subroutine sd_verlet_ciccotti_setup
2301 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's 
2302 c version 
2303       implicit real*8 (a-h,o-z)
2304       include 'DIMENSIONS'
2305 #ifdef MPI
2306       include 'mpif.h'
2307 #endif
2308       include 'COMMON.CONTROL'
2309       include 'COMMON.VAR'
2310       include 'COMMON.MD'
2311 #ifndef LANG0
2312       include 'COMMON.LANGEVIN'
2313 #else
2314       include 'COMMON.LANGEVIN.lang0'
2315 #endif
2316       include 'COMMON.CHAIN'
2317       include 'COMMON.DERIV'
2318       include 'COMMON.GEO'
2319       include 'COMMON.LOCAL'
2320       include 'COMMON.INTERACT'
2321       include 'COMMON.IOUNITS'
2322       include 'COMMON.NAMES'
2323       include 'COMMON.TIME1'
2324       double precision emgdt(MAXRES6),
2325      & pterm,vterm,rho,rhoc,vsig,
2326      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2327      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2328      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2329       logical lprn /.false./
2330       double precision zero /1.0d-8/, gdt_radius /0.05d0/ 
2331       double precision ktm
2332 #ifdef MPI
2333       tt0 = MPI_Wtime()
2334 #else
2335       tt0 = tcpu()
2336 #endif
2337 c
2338 c AL 8/17/04 Code adapted from tinker
2339 c
2340 c Get the frictional and random terms for stochastic dynamics in the
2341 c eigenspace of mass-scaled UNRES friction matrix
2342 c
2343       do i = 1, dimen
2344             write (iout,*) "i",i," fricgam",fricgam(i)
2345             gdt = fricgam(i) * d_time
2346 c
2347 c Stochastic dynamics reduces to simple MD for zero friction
2348 c
2349             if (gdt .le. zero) then
2350                pfric_vec(i) = 1.0d0
2351                vfric_vec(i) = d_time
2352                afric_vec(i) = 0.5d0*d_time*d_time
2353                prand_vec(i) = afric_vec(i)
2354                vrand_vec2(i) = vfric_vec(i)
2355 c
2356 c Analytical expressions when friction coefficient is large
2357 c
2358             else 
2359                egdt = dexp(-gdt)
2360                pfric_vec(i) = egdt
2361                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2362                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2363                prand_vec(i) = afric_vec(i)
2364                vrand_vec2(i) = vfric_vec(i)
2365 c
2366 c Compute the scaling factors of random terms for the nonzero friction case
2367 c
2368 c               ktm = 0.5d0*d_time/fricgam(i)
2369 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2370 c               vsig = dsqrt(ktm*vterm)
2371 c               prand_vec(i) = psig*afric_vec(i) 
2372 c               vrand_vec2(i) = vsig*vfric_vec(i)
2373             end if
2374       end do
2375       if (lprn) then
2376       write (iout,*) 
2377      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2378      &  " vrand_vec2"
2379       do i=1,dimen
2380         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2381      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2382       enddo
2383       endif
2384 c
2385 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2386 c
2387       call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
2388       call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
2389       call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
2390       call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
2391       call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
2392 #ifdef MPI
2393       t_sdsetup=t_sdsetup+MPI_Wtime()
2394 #else
2395       t_sdsetup=t_sdsetup+tcpu()-tt0
2396 #endif
2397       return
2398       end
2399 c-------------------------------------------------------------      
2400       subroutine sd_verlet1_ciccotti
2401 c Applying stochastic velocity Verlet algorithm - step 1 to velocities        
2402       implicit real*8 (a-h,o-z)
2403       include 'DIMENSIONS'
2404 #ifdef MPI
2405       include 'mpif.h'
2406 #endif
2407       include 'COMMON.CONTROL'
2408       include 'COMMON.VAR'
2409       include 'COMMON.MD'
2410 #ifndef LANG0
2411       include 'COMMON.LANGEVIN'
2412 #else
2413       include 'COMMON.LANGEVIN.lang0'
2414 #endif
2415       include 'COMMON.CHAIN'
2416       include 'COMMON.DERIV'
2417       include 'COMMON.GEO'
2418       include 'COMMON.LOCAL'
2419       include 'COMMON.INTERACT'
2420       include 'COMMON.IOUNITS'
2421       include 'COMMON.NAMES'
2422       double precision stochforcvec(MAXRES6)
2423       common /stochcalc/ stochforcvec
2424       logical lprn /.false./
2425
2426 c      write (iout,*) "dc_old"
2427 c      do i=0,nres
2428 c        write (iout,'(i5,3f10.5,5x,3f10.5)') 
2429 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2430 c      enddo
2431       do j=1,3
2432         dc_work(j)=dc_old(j,0)
2433         d_t_work(j)=d_t_old(j,0)
2434         d_a_work(j)=d_a_old(j,0)
2435       enddo
2436       ind=3
2437       do i=nnt,nct-1
2438         do j=1,3
2439           dc_work(ind+j)=dc_old(j,i)
2440           d_t_work(ind+j)=d_t_old(j,i)
2441           d_a_work(ind+j)=d_a_old(j,i)
2442         enddo
2443         ind=ind+3
2444       enddo
2445       do i=nnt,nct
2446         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2447           do j=1,3
2448             dc_work(ind+j)=dc_old(j,i+nres)
2449             d_t_work(ind+j)=d_t_old(j,i+nres)
2450             d_a_work(ind+j)=d_a_old(j,i+nres)
2451           enddo
2452           ind=ind+3
2453         endif
2454       enddo
2455
2456 #ifndef LANG0
2457       if (lprn) then
2458       write (iout,*) 
2459      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2460      &  " vrand_mat2"
2461       do i=1,dimen
2462         do j=1,dimen
2463           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2464      &      vfric_mat(i,j),afric_mat(i,j),
2465      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2466         enddo
2467       enddo
2468       endif
2469       do i=1,dimen
2470         ddt1=0.0d0
2471         ddt2=0.0d0
2472         do j=1,dimen
2473           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2474      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2475           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2476           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2477         enddo
2478         d_t_work_new(i)=ddt1+0.5d0*ddt2
2479         d_t_work(i)=ddt1+ddt2
2480       enddo
2481 #endif
2482       do j=1,3
2483         dc(j,0)=dc_work(j)
2484         d_t(j,0)=d_t_work(j)
2485       enddo
2486       ind=3     
2487       do i=nnt,nct-1    
2488         do j=1,3
2489           dc(j,i)=dc_work(ind+j)
2490           d_t(j,i)=d_t_work(ind+j)
2491         enddo
2492         ind=ind+3
2493       enddo
2494       do i=nnt,nct
2495         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2496           inres=i+nres
2497           do j=1,3
2498             dc(j,inres)=dc_work(ind+j)
2499             d_t(j,inres)=d_t_work(ind+j)
2500           enddo
2501           ind=ind+3
2502         endif      
2503       enddo 
2504       return
2505       end
2506 c--------------------------------------------------------------------------
2507       subroutine sd_verlet2_ciccotti
2508 c  Calculating the adjusted velocities for accelerations
2509       implicit real*8 (a-h,o-z)
2510       include 'DIMENSIONS'
2511       include 'COMMON.CONTROL'
2512       include 'COMMON.VAR'
2513       include 'COMMON.MD'
2514 #ifndef LANG0
2515       include 'COMMON.LANGEVIN'
2516 #else
2517       include 'COMMON.LANGEVIN.lang0'
2518 #endif
2519       include 'COMMON.CHAIN'
2520       include 'COMMON.DERIV'
2521       include 'COMMON.GEO'
2522       include 'COMMON.LOCAL'
2523       include 'COMMON.INTERACT'
2524       include 'COMMON.IOUNITS'
2525       include 'COMMON.NAMES'
2526       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2527       common /stochcalc/ stochforcvec
2528 c
2529 c Compute the stochastic forces which contribute to velocity change
2530 c
2531       call stochastic_force(stochforcvecV)
2532 #ifndef LANG0
2533       do i=1,dimen
2534         ddt1=0.0d0
2535         ddt2=0.0d0
2536         do j=1,dimen
2537
2538           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2539 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2540           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2541         enddo
2542         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2543       enddo
2544 #endif
2545       do j=1,3
2546         d_t(j,0)=d_t_work(j)
2547       enddo
2548       ind=3
2549       do i=nnt,nct-1
2550         do j=1,3
2551           d_t(j,i)=d_t_work(ind+j)
2552         enddo
2553         ind=ind+3
2554       enddo
2555       do i=nnt,nct
2556         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2557           inres=i+nres
2558           do j=1,3
2559             d_t(j,inres)=d_t_work(ind+j)
2560           enddo
2561           ind=ind+3
2562         endif
2563       enddo 
2564       return
2565       end