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