added source code
[unres.git] / source / unres / src_MD / src / md-diff / mts / sd_verlet1_ciccotti.f
1 c-------------------------------------------------------------      
2       subroutine sd_verlet1_ciccotti
3 c Applying stochastic velocity Verlet algorithm - step 1 to velocities        
4       implicit real*8 (a-h,o-z)
5       include 'DIMENSIONS'
6 #ifdef MPI
7       include 'mpif.h'
8 #endif
9       include 'COMMON.CONTROL'
10       include 'COMMON.VAR'
11       include 'COMMON.MD'
12 #ifndef LANG0
13       include 'COMMON.LANGEVIN'
14 #else
15       include 'COMMON.LANGEVIN.lang0'
16 #endif
17       include 'COMMON.CHAIN'
18       include 'COMMON.DERIV'
19       include 'COMMON.GEO'
20       include 'COMMON.LOCAL'
21       include 'COMMON.INTERACT'
22       include 'COMMON.IOUNITS'
23       include 'COMMON.NAMES'
24       double precision stochforcvec(MAXRES6)
25       common /stochcalc/ stochforcvec
26       logical lprn /.false./
27
28 c      write (iout,*) "dc_old"
29 c      do i=0,nres
30 c        write (iout,'(i5,3f10.5,5x,3f10.5)') 
31 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
32 c      enddo
33       do j=1,3
34         dc_work(j)=dc_old(j,0)
35         d_t_work(j)=d_t_old(j,0)
36         d_a_work(j)=d_a_old(j,0)
37       enddo
38       ind=3
39       do i=nnt,nct-1
40         do j=1,3
41           dc_work(ind+j)=dc_old(j,i)
42           d_t_work(ind+j)=d_t_old(j,i)
43           d_a_work(ind+j)=d_a_old(j,i)
44         enddo
45         ind=ind+3
46       enddo
47       do i=nnt,nct
48         if (itype(i).ne.10) then
49           do j=1,3
50             dc_work(ind+j)=dc_old(j,i+nres)
51             d_t_work(ind+j)=d_t_old(j,i+nres)
52             d_a_work(ind+j)=d_a_old(j,i+nres)
53           enddo
54           ind=ind+3
55         endif
56       enddo
57
58 #ifndef LANG0
59       if (lprn) then
60       write (iout,*) 
61      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
62      &  " vrand_mat2"
63       do i=1,dimen
64         do j=1,dimen
65           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
66      &      vfric_mat(i,j),afric_mat(i,j),
67      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
68         enddo
69       enddo
70       endif
71       do i=1,dimen
72         ddt1=0.0d0
73         ddt2=0.0d0
74         do j=1,dimen
75           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
76      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
77           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
78           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
79         enddo
80         d_t_work_new(i)=ddt1+0.5d0*ddt2
81         d_t_work(i)=ddt1+ddt2
82       enddo
83 #endif
84       do j=1,3
85         dc(j,0)=dc_work(j)
86         d_t(j,0)=d_t_work(j)
87       enddo
88       ind=3     
89       do i=nnt,nct-1    
90         do j=1,3
91           dc(j,i)=dc_work(ind+j)
92           d_t(j,i)=d_t_work(ind+j)
93         enddo
94         ind=ind+3
95       enddo
96       do i=nnt,nct
97         if (itype(i).ne.10) then
98           inres=i+nres
99           do j=1,3
100             dc(j,inres)=dc_work(ind+j)
101             d_t(j,inres)=d_t_work(ind+j)
102           enddo
103           ind=ind+3
104         endif      
105       enddo 
106       return
107       end