Merge branch 'bartek' of mmka.chem.univ.gda.pl:unres into bartek
[unres.git] / source / unres / src_MD / moments.f
1       subroutine inertia_tensor
2 c Calculating the intertia tensor for the entire protein in order to
3 c remove the perpendicular components of velocity matrix which cause
4 c the molecule to rotate.       
5        implicit real*8 (a-h,o-z)
6        include 'DIMENSIONS'
7        include 'COMMON.CONTROL'
8        include 'COMMON.VAR'
9        include 'COMMON.MD'
10        include 'COMMON.CHAIN'
11        include 'COMMON.DERIV'
12        include 'COMMON.GEO'
13        include 'COMMON.LOCAL'
14        include 'COMMON.INTERACT'
15        include 'COMMON.IOUNITS'
16        include 'COMMON.NAMES'
17       
18       double precision Im(3,3),Imcp(3,3),cm(3),pr(3),M_SC,
19      & eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),
20      & vpp(3,0:MAXRES),vs_p(3),pr1(3,3),
21      & pr2(3,3),pp(3),incr(3),v(3),mag,mag2 
22       common /gucio/ cm
23       integer iti,inres 
24         do i=1,3
25            do j=1,3
26              Im(i,j)=0.0d0
27              pr1(i,j)=0.0d0
28              pr2(i,j)=0.0d0                  
29            enddo
30            L(i)=0.0d0
31            cm(i)=0.0d0
32            vrot(i)=0.0d0                   
33         enddo
34 c   calculating the center of the mass of the protein                                   
35         do i=nnt,nct-1
36           do j=1,3
37             cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
38           enddo
39         enddo
40         do j=1,3
41          cm(j)=mp*cm(j)
42         enddo
43         M_SC=0.0d0                              
44         do i=nnt,nct
45            iti=itype(i)          
46            M_SC=M_SC+msc(iti)
47            inres=i+nres
48            do j=1,3
49             cm(j)=cm(j)+msc(iti)*c(j,inres)         
50            enddo
51         enddo
52         do j=1,3
53           cm(j)=cm(j)/(M_SC+(nct-nnt)*mp)
54         enddo
55        
56         do i=nnt,nct-1
57           do j=1,3
58             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
59           enddo
60           Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3))
61           Im(1,2)=Im(1,2)-mp*pr(1)*pr(2)
62           Im(1,3)=Im(1,3)-mp*pr(1)*pr(3)
63           Im(2,3)=Im(2,3)-mp*pr(2)*pr(3)        
64           Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1))
65           Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2))
66         enddo                   
67         
68         do i=nnt,nct    
69            iti=itype(i)
70            inres=i+nres
71            do j=1,3
72              pr(j)=c(j,inres)-cm(j)         
73            enddo
74           Im(1,1)=Im(1,1)+msc(iti)*(pr(2)*pr(2)+pr(3)*pr(3))
75           Im(1,2)=Im(1,2)-msc(iti)*pr(1)*pr(2)
76           Im(1,3)=Im(1,3)-msc(iti)*pr(1)*pr(3)
77           Im(2,3)=Im(2,3)-msc(iti)*pr(2)*pr(3)  
78           Im(2,2)=Im(2,2)+msc(iti)*(pr(3)*pr(3)+pr(1)*pr(1))
79           Im(3,3)=Im(3,3)+msc(iti)*(pr(1)*pr(1)+pr(2)*pr(2))               
80         enddo
81           
82         do i=nnt,nct-1
83           Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))*       
84      &    vbld(i+1)*vbld(i+1)*0.25d0
85           Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))*
86      &    vbld(i+1)*vbld(i+1)*0.25d0              
87           Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))*
88      &    vbld(i+1)*vbld(i+1)*0.25d0      
89           Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))*
90      &    vbld(i+1)*vbld(i+1)*0.25d0            
91           Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))*
92      &    vbld(i+1)*vbld(i+1)*0.25d0      
93           Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))*
94      &    vbld(i+1)*vbld(i+1)*0.25d0
95         enddo
96         
97                                 
98         do i=nnt,nct
99          if (itype(i).ne.10) then
100            iti=itype(i)          
101            inres=i+nres
102           Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)*
103      &    dc_norm(1,inres))*vbld(inres)*vbld(inres)
104           Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)*
105      &    dc_norm(2,inres))*vbld(inres)*vbld(inres)
106           Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)*
107      &    dc_norm(3,inres))*vbld(inres)*vbld(inres)
108           Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)*
109      &    dc_norm(3,inres))*vbld(inres)*vbld(inres)     
110           Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)*
111      &    dc_norm(2,inres))*vbld(inres)*vbld(inres)
112           Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)*
113      &    dc_norm(3,inres))*vbld(inres)*vbld(inres)
114          endif
115         enddo
116         
117         call angmom(cm,L)
118 c        write(iout,*) "The angular momentum before adjustment:"
119 c        write(iout,*) (L(j),j=1,3)                                                                                                                                                                                                                     
120         
121         Im(2,1)=Im(1,2)
122         Im(3,1)=Im(1,3)
123         Im(3,2)=Im(2,3)
124       
125 c  Copying the Im matrix for the djacob subroutine
126         do i=1,3
127           do j=1,3
128             Imcp(i,j)=Im(i,j)
129             Id(i,j)=0.0d0           
130           enddo
131         enddo
132                                                               
133 c   Finding the eigenvectors and eignvalues of the inertia tensor
134        call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
135 c       write (iout,*) "Eigenvalues & Eigenvectors"
136 c       write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
137 c       write (iout,*)
138 c       do i=1,3
139 c         write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
140 c       enddo
141 c   Constructing the diagonalized matrix
142        do i=1,3
143          if (dabs(eigval(i)).gt.1.0d-15) then
144            Id(i,i)=1.0d0/eigval(i)
145          else
146            Id(i,i)=0.0d0
147          endif
148        enddo
149         do i=1,3
150            do j=1,3
151               Imcp(i,j)=eigvec(j,i)
152            enddo
153         enddo    
154         do i=1,3
155            do j=1,3
156               do k=1,3   
157                  pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
158               enddo
159            enddo
160         enddo
161         do i=1,3
162            do j=1,3
163               do k=1,3   
164                  pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
165               enddo
166            enddo
167         enddo
168 c  Calculating the total rotational velocity of the molecule
169        do i=1,3    
170          do j=1,3
171            vrot(i)=vrot(i)+pr2(i,j)*L(j)
172          enddo
173        enddo    
174 c   Resetting the velocities
175        do i=nnt,nct-1
176          call vecpr(vrot(1),dc(1,i),vp)  
177          do j=1,3
178            d_t(j,i)=d_t(j,i)-vp(j)
179           enddo
180         enddo
181         do i=nnt,nct 
182         if(itype(i).ne.10) then
183            inres=i+nres
184            call vecpr(vrot(1),dc(1,inres),vp)                    
185            do j=1,3
186              d_t(j,inres)=d_t(j,inres)-vp(j)
187            enddo
188         endif
189        enddo
190        call angmom(cm,L)
191 c       write(iout,*) "The angular momentum after adjustment:"
192 c       write(iout,*) (L(j),j=1,3)                                                                                                                                                                                                                      
193        return
194        end 
195 c----------------------------------------------------------------------------
196        subroutine angmom(cm,L)
197        implicit real*8 (a-h,o-z)
198        include 'DIMENSIONS'
199        include 'COMMON.CONTROL'
200        include 'COMMON.VAR'
201        include 'COMMON.MD'
202        include 'COMMON.CHAIN'
203        include 'COMMON.DERIV'
204        include 'COMMON.GEO'
205        include 'COMMON.LOCAL'
206        include 'COMMON.INTERACT'
207        include 'COMMON.IOUNITS'
208        include 'COMMON.NAMES'
209       
210       double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),
211      &  pp(3)
212       integer iti,inres 
213 c  Calculate the angular momentum
214        do j=1,3
215           L(j)=0.0d0
216        enddo
217        do j=1,3
218           incr(j)=d_t(j,0)
219        enddo                   
220        do i=nnt,nct-1
221           do j=1,3
222             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
223           enddo
224           do j=1,3
225             v(j)=incr(j)+0.5d0*d_t(j,i)
226           enddo
227           do j=1,3
228             incr(j)=incr(j)+d_t(j,i)
229           enddo         
230           call vecpr(pr(1),v(1),vp)
231           do j=1,3
232             L(j)=L(j)+mp*vp(j)
233           enddo
234           do j=1,3
235              pr(j)=0.5d0*dc(j,i)
236              pp(j)=0.5d0*d_t(j,i)                 
237           enddo
238          call vecpr(pr(1),pp(1),vp)
239          do j=1,3                
240              L(j)=L(j)+Ip*vp(j)  
241           enddo
242         enddo
243         do j=1,3
244           incr(j)=d_t(j,0)
245         enddo   
246         do i=nnt,nct
247          iti=itype(i)    
248          inres=i+nres
249          do j=1,3
250            pr(j)=c(j,inres)-cm(j)           
251          enddo
252          if (itype(i).ne.10) then
253            do j=1,3
254              v(j)=incr(j)+d_t(j,inres)
255            enddo
256          else
257            do j=1,3
258              v(j)=incr(j)
259            enddo
260          endif
261          call vecpr(pr(1),v(1),vp)
262 c         write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
263 c     &     " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
264          do j=1,3
265             L(j)=L(j)+msc(iti)*vp(j)
266          enddo
267 c         write (iout,*) "L",(l(j),j=1,3)
268          if (itype(i).ne.10) then
269            do j=1,3
270             v(j)=incr(j)+d_t(j,inres)
271            enddo
272            call vecpr(dc(1,inres),d_t(1,inres),vp)
273            do j=1,3                        
274              L(j)=L(j)+Isc(iti)*vp(j)    
275           enddo                    
276          endif
277          do j=1,3
278              incr(j)=incr(j)+d_t(j,i)
279          enddo
280        enddo
281       return
282       end
283 c------------------------------------------------------------------------------
284        subroutine vcm_vel(vcm)
285        implicit real*8 (a-h,o-z)
286        include 'DIMENSIONS'
287        include 'COMMON.VAR'
288        include 'COMMON.MD'
289        include 'COMMON.CHAIN'
290        include 'COMMON.DERIV'
291        include 'COMMON.GEO'
292        include 'COMMON.LOCAL'
293        include 'COMMON.INTERACT'
294        include 'COMMON.IOUNITS'
295        double precision vcm(3),vv(3),summas,amas
296        do j=1,3
297          vcm(j)=0.0d0
298          vv(j)=d_t(j,0)
299        enddo
300        summas=0.0d0
301        do i=nnt,nct
302          if (i.lt.nct) then
303            summas=summas+mp
304            do j=1,3
305              vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
306            enddo
307          endif
308          amas=msc(itype(i))
309          summas=summas+amas              
310          if (itype(i).ne.10) then
311            do j=1,3
312              vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
313            enddo
314          else
315            do j=1,3
316              vcm(j)=vcm(j)+amas*vv(j)
317            enddo
318          endif
319          do j=1,3
320            vv(j)=vv(j)+d_t(j,i)
321          enddo
322        enddo 
323 c       write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
324        do j=1,3
325          vcm(j)=vcm(j)/summas
326        enddo
327        return
328        end