Merge branch 'homology' of mmka.chem.univ.gda.pl:unres into homology
[unres.git] / source / unres / src_MD-M / q_measure.F
1       double precision function qwolynes(seg1,seg2,flag,seg3,seg4)
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'COMMON.IOUNITS'
5       include 'COMMON.CHAIN' 
6       include 'COMMON.INTERACT'
7       include 'COMMON.VAR'
8       integer i,j,jl,k,l,il,kl,nl,np,ip,kp,seg1,seg2,seg3,seg4,
9      & secseg
10       integer nsep /3/
11       double precision dist,qm
12       double precision qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
13       logical lprn /.false./
14       logical flag
15       double precision sigm,x
16       sigm(x)=0.25d0*x
17       do kkk=1,nperm
18       qq = 0.0d0
19       nl=0 
20        if(flag) then
21         do il=seg1+nsep,seg2
22           do jl=seg1,il-nsep
23             nl=nl+1
24             d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+
25      &                 (cref(2,jl,kkk)-cref(2,il,kkk))**2+
26      &                 (cref(3,jl,kkk)-cref(3,il,kkk))**2)
27             dij=dist(il,jl)
28             qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
29             if (itype(il).ne.10 .or. itype(jl).ne.10) then
30               nl=nl+1
31               d0ijCM=dsqrt(
32      &               (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+
33      &               (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+
34      &               (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
35               dijCM=dist(il+nres,jl+nres)
36               qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
37             endif
38             qq = qq+qqij+qqijCM
39           enddo
40         enddo   
41         qq = qq/nl
42       else
43       do il=seg1,seg2
44         if((seg3-il).lt.3) then
45              secseg=il+3
46         else
47              secseg=seg3
48         endif 
49           do jl=secseg,seg4
50             nl=nl+1
51             d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+
52      &                 (cref(2,jl,kkk)-cref(2,il,kkk))**2+
53      &                 (cref(3,jl,kkk)-cref(3,il,kkk))**2)
54             dij=dist(il,jl)
55             qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
56             if (itype(il).ne.10 .or. itype(jl).ne.10) then
57               nl=nl+1
58               d0ijCM=dsqrt(
59      &               (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+
60      &               (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+
61      &               (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
62               dijCM=dist(il+nres,jl+nres)
63               qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
64             endif
65             qq = qq+qqij+qqijCM
66           enddo
67         enddo
68       qq = qq/nl
69       endif
70       if (qqmax.le.qq) qqmax=qq
71       enddo
72       qwolynes=1.0d0-qqmax
73       return 
74       end
75 c-------------------------------------------------------------------
76       subroutine qwolynes_prim(seg1,seg2,flag,seg3,seg4)
77       implicit real*8 (a-h,o-z)
78       include 'DIMENSIONS'
79       include 'COMMON.IOUNITS'
80       include 'COMMON.CHAIN' 
81       include 'COMMON.INTERACT'
82       include 'COMMON.VAR'
83       include 'COMMON.MD'
84       integer i,j,jl,k,l,il,nl,seg1,seg2,seg3,seg4,
85      & secseg
86       integer nsep /3/
87       double precision dist
88       double precision dij,d0ij,dijCM,d0ijCM
89       logical lprn /.false./
90       logical flag
91       double precision sigm,x,sim,dd0,fac,ddqij
92       sigm(x)=0.25d0*x
93       do kkk=1,nperm 
94       do i=0,nres
95         do j=1,3
96           dqwol(j,i)=0.0d0
97           dxqwol(j,i)=0.0d0       
98         enddo
99       enddo
100       nl=0 
101        if(flag) then
102         do il=seg1+nsep,seg2
103           do jl=seg1,il-nsep
104             nl=nl+1
105             d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+
106      &                 (cref(2,jl,kkk)-cref(2,il,kkk))**2+
107      &                 (cref(3,jl,kkk)-cref(3,il,kkk))**2)
108             dij=dist(il,jl)
109             sim = 1.0d0/sigm(d0ij)
110             sim = sim*sim
111             dd0 = dij-d0ij
112             fac = dd0*sim/dij*dexp(-0.5d0*dd0*dd0*sim)
113             do k=1,3
114               ddqij = (c(k,il)-c(k,jl))*fac
115               dqwol(k,il)=dqwol(k,il)+ddqij
116               dqwol(k,jl)=dqwol(k,jl)-ddqij
117             enddo
118                      
119             if (itype(il).ne.10 .or. itype(jl).ne.10) then
120               nl=nl+1
121               d0ijCM=dsqrt(
122      &               (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+
123      &               (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+
124      &               (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
125               dijCM=dist(il+nres,jl+nres)
126               sim = 1.0d0/sigm(d0ijCM)
127               sim = sim*sim
128               dd0=dijCM-d0ijCM
129               fac=dd0*sim/dijCM*dexp(-0.5d0*dd0*dd0*sim)
130               do k=1,3
131                 ddqij = (c(k,il+nres)-c(k,jl+nres))*fac
132                 dxqwol(k,il)=dxqwol(k,il)+ddqij
133                 dxqwol(k,jl)=dxqwol(k,jl)-ddqij
134               enddo
135             endif           
136           enddo
137         enddo   
138        else
139         do il=seg1,seg2
140         if((seg3-il).lt.3) then
141              secseg=il+3
142         else
143              secseg=seg3
144         endif 
145           do jl=secseg,seg4
146             nl=nl+1
147             d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+
148      &                 (cref(2,jl,kkk)-cref(2,il,kkk))**2+
149      &                 (cref(3,jl,kkk)-cref(3,il,kkk))**2)
150             dij=dist(il,jl)
151             sim = 1.0d0/sigm(d0ij)
152             sim = sim*sim
153             dd0 = dij-d0ij
154             fac = dd0*sim/dij*dexp(-0.5d0*dd0*dd0*sim)
155             do k=1,3
156               ddqij = (c(k,il)-c(k,jl))*fac
157               dqwol(k,il)=dqwol(k,il)+ddqij
158               dqwol(k,jl)=dqwol(k,jl)-ddqij
159             enddo
160             if (itype(il).ne.10 .or. itype(jl).ne.10) then
161               nl=nl+1
162               d0ijCM=dsqrt(
163      &               (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+
164      &               (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+
165      &               (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
166               dijCM=dist(il+nres,jl+nres)
167               sim = 1.0d0/sigm(d0ijCM)
168               sim=sim*sim
169               dd0 = dijCM-d0ijCM
170               fac = dd0*sim/dijCM*dexp(-0.5d0*dd0*dd0*sim)
171               do k=1,3
172                ddqij = (c(k,il+nres)-c(k,jl+nres))*fac             
173                dxqwol(k,il)=dxqwol(k,il)+ddqij
174                dxqwol(k,jl)=dxqwol(k,jl)-ddqij  
175               enddo
176             endif 
177           enddo
178         enddo                
179       endif
180       enddo
181        do i=0,nres
182          do j=1,3
183            dqwol(j,i)=dqwol(j,i)/nl
184            dxqwol(j,i)=dxqwol(j,i)/nl
185          enddo
186        enddo                                                             
187       return 
188       end
189 c-------------------------------------------------------------------
190       subroutine qwol_num(seg1,seg2,flag,seg3,seg4)
191       implicit real*8 (a-h,o-z)
192       include 'DIMENSIONS'
193       include 'COMMON.IOUNITS'
194       include 'COMMON.CHAIN' 
195       include 'COMMON.INTERACT'
196       include 'COMMON.VAR'
197       integer seg1,seg2,seg3,seg4
198       logical flag
199       double precision qwolan(3,0:maxres),cdummy(3,0:maxres2),
200      & qwolxan(3,0:maxres),q1,q2
201       double precision delta /1.0d-10/
202       do i=0,nres
203         do j=1,3
204           q1=qwolynes(seg1,seg2,flag,seg3,seg4)
205           cdummy(j,i)=c(j,i)
206           c(j,i)=c(j,i)+delta
207           q2=qwolynes(seg1,seg2,flag,seg3,seg4)
208           qwolan(j,i)=(q2-q1)/delta
209           c(j,i)=cdummy(j,i)
210         enddo
211       enddo
212       do i=0,nres
213         do j=1,3
214           q1=qwolynes(seg1,seg2,flag,seg3,seg4)
215           cdummy(j,i+nres)=c(j,i+nres)
216           c(j,i+nres)=c(j,i+nres)+delta
217           q2=qwolynes(seg1,seg2,flag,seg3,seg4)
218           qwolxan(j,i)=(q2-q1)/delta
219           c(j,i+nres)=cdummy(j,i+nres)
220         enddo
221       enddo  
222 c      write(iout,*) "Numerical Q carteisan gradients backbone: "
223 c      do i=0,nct
224 c        write(iout,'(i5,3e15.5)') i, (qwolan(j,i),j=1,3)
225 c      enddo
226 c      write(iout,*) "Numerical Q carteisan gradients side-chain: "
227 c      do i=0,nct
228 c        write(iout,'(i5,3e15.5)') i, (qwolxan(j,i),j=1,3)
229 c      enddo 
230       return
231       end
232 c------------------------------------------------------------------------  
233       subroutine EconstrQ
234 c     MD with umbrella_sampling using Wolyne's distance measure as a constraint
235       implicit real*8 (a-h,o-z)
236       include 'DIMENSIONS'
237       include 'COMMON.CONTROL'
238       include 'COMMON.VAR'
239       include 'COMMON.MD'
240 #ifndef LANG0
241       include 'COMMON.LANGEVIN'
242 #else
243       include 'COMMON.LANGEVIN.lang0'
244 #endif
245       include 'COMMON.CHAIN'
246       include 'COMMON.DERIV'
247       include 'COMMON.GEO'
248       include 'COMMON.LOCAL'
249       include 'COMMON.INTERACT'
250       include 'COMMON.IOUNITS'
251       include 'COMMON.NAMES'
252       include 'COMMON.TIME1'
253       double precision uzap1,uzap2,hm1,hm2,hmnum
254       double precision ucdelan,dUcartan(3,0:MAXRES)
255      & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
256      &  duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
257       integer kstart,kend,lstart,lend,idummy
258       double precision delta /1.0d-7/
259       do i=0,nres
260          do j=1,3
261             duconst(j,i)=0.0d0
262             dudconst(j,i)=0.0d0     
263             duxconst(j,i)=0.0d0
264             dudxconst(j,i)=0.0d0             
265          enddo
266       enddo
267       Uconst=0.0d0
268       do i=1,nfrag
269          qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
270      &    ,idummy,idummy)
271          Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
272 c Calculating the derivatives of Constraint energy with respect to Q
273          Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),
274      &     qinfrag(i,iset))
275 c         hm1=harmonic(qfrag(i,iset),qinfrag(i,iset))
276 c        hm2=harmonic(qfrag(i,iset)+delta,qinfrag(i,iset))
277 c         hmnum=(hm2-hm1)/delta          
278 c         write(iout,*) "harmonicprim frag",harmonicprim(qfrag(i,iset),
279 c     &   qinfrag(i,iset))
280 c         write(iout,*) "harmonicnum frag", hmnum                
281 c Calculating the derivatives of Q with respect to cartesian coordinates
282          call qwolynes_prim(ifrag(1,i,iset),ifrag(2,i,iset),.true.
283      &   ,idummy,idummy)
284 c         write(iout,*) "dqwol "
285 c         do ii=1,nres
286 c          write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
287 c         enddo
288 c         write(iout,*) "dxqwol "
289 c         do ii=1,nres
290 c           write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
291 c         enddo
292 c Calculating numerical gradients of dU/dQi and dQi/dxi
293 c        call qwol_num(ifrag(1,i,iset),ifrag(2,i,iset),.true.
294 c     &  ,idummy,idummy)
295 c  The gradients of Uconst in Cs
296          do ii=0,nres
297             do j=1,3
298                duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
299                dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
300             enddo
301          enddo
302       enddo     
303       do i=1,npair
304          kstart=ifrag(1,ipair(1,i,iset),iset)
305          kend=ifrag(2,ipair(1,i,iset),iset)
306          lstart=ifrag(1,ipair(2,i,iset),iset)
307          lend=ifrag(2,ipair(2,i,iset),iset)
308          qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
309          Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
310 c  Calculating dU/dQ
311          Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
312 c         hm1=harmonic(qpair(i),qinpair(i,iset))
313 c        hm2=harmonic(qpair(i)+delta,qinpair(i,iset))
314 c         hmnum=(hm2-hm1)/delta          
315 c         write(iout,*) "harmonicprim pair ",harmonicprim(qpair(i),
316 c     &   qinpair(i,iset))
317 c         write(iout,*) "harmonicnum pair ", hmnum       
318 c Calculating dQ/dXi
319          call qwolynes_prim(kstart,kend,.false.
320      &   ,lstart,lend)
321 c         write(iout,*) "dqwol "
322 c         do ii=1,nres
323 c          write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
324 c         enddo
325 c         write(iout,*) "dxqwol "
326 c         do ii=1,nres
327 c          write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
328 c        enddo
329 c Calculating numerical gradients
330 c        call qwol_num(kstart,kend,.false.
331 c     &  ,lstart,lend)
332 c The gradients of Uconst in Cs
333          do ii=0,nres
334             do j=1,3
335                duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
336                dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
337             enddo
338          enddo
339       enddo
340 c      write(iout,*) "Uconst inside subroutine ", Uconst
341 c Transforming the gradients from Cs to dCs for the backbone
342       do i=0,nres
343          do j=i+1,nres
344            do k=1,3
345              dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
346            enddo
347          enddo
348       enddo
349 c  Transforming the gradients from Cs to dCs for the side chains      
350       do i=1,nres
351          do j=1,3
352            dudxconst(j,i)=duxconst(j,i)
353          enddo
354       enddo                      
355 c      write(iout,*) "dU/ddc backbone "
356 c       do ii=0,nres
357 c        write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
358 c      enddo      
359 c      write(iout,*) "dU/ddX side chain "
360 c      do ii=1,nres
361 c            write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
362 c      enddo
363 c Calculating numerical gradients of dUconst/ddc and dUconst/ddx
364 c      call dEconstrQ_num      
365       return
366       end
367 c-----------------------------------------------------------------------
368       subroutine dEconstrQ_num
369 c Calculating numerical dUconst/ddc and dUconst/ddx      
370       implicit real*8 (a-h,o-z)
371       include 'DIMENSIONS'
372       include 'COMMON.CONTROL'
373       include 'COMMON.VAR'
374       include 'COMMON.MD'
375 #ifndef LANG0
376       include 'COMMON.LANGEVIN'
377 #else
378       include 'COMMON.LANGEVIN.lang0'
379 #endif
380       include 'COMMON.CHAIN'
381       include 'COMMON.DERIV'
382       include 'COMMON.GEO'
383       include 'COMMON.LOCAL'
384       include 'COMMON.INTERACT'
385       include 'COMMON.IOUNITS'
386       include 'COMMON.NAMES'
387       include 'COMMON.TIME1'
388       double precision uzap1,uzap2
389       double precision dUcartan(3,0:MAXRES)
390      & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES)
391       integer kstart,kend,lstart,lend,idummy
392       double precision delta /1.0d-7/
393 c     For the backbone
394       do i=0,nres-1
395          do j=1,3
396             dUcartan(j,i)=0.0d0
397             cdummy(j,i)=dc(j,i)
398             dc(j,i)=dc(j,i)+delta
399             call chainbuild_cart
400             uzap2=0.0d0
401             do ii=1,nfrag
402              qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
403      &         ,idummy,idummy)
404                uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
405      &          qinfrag(ii,iset))
406             enddo
407             do ii=1,npair
408                kstart=ifrag(1,ipair(1,ii,iset),iset)
409                kend=ifrag(2,ipair(1,ii,iset),iset)
410                lstart=ifrag(1,ipair(2,ii,iset),iset)
411                lend=ifrag(2,ipair(2,ii,iset),iset)
412                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
413                uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),
414      &           qinpair(ii,iset))
415             enddo
416             dc(j,i)=cdummy(j,i)
417             call chainbuild_cart
418             uzap1=0.0d0
419              do ii=1,nfrag
420              qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
421      &         ,idummy,idummy)
422                uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),
423      &          qinfrag(ii,iset))
424             enddo
425             do ii=1,npair
426                kstart=ifrag(1,ipair(1,ii,iset),iset)
427                kend=ifrag(2,ipair(1,ii,iset),iset)
428                lstart=ifrag(1,ipair(2,ii,iset),iset)
429                lend=ifrag(2,ipair(2,ii,iset),iset)
430                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
431                uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),
432      &          qinpair(ii,iset))
433             enddo
434             ducartan(j,i)=(uzap2-uzap1)/(delta)     
435          enddo
436       enddo
437 c Calculating numerical gradients for dU/ddx
438       do i=0,nres-1
439          duxcartan(j,i)=0.0d0
440          do j=1,3
441             cdummy(j,i)=dc(j,i+nres)
442             dc(j,i+nres)=dc(j,i+nres)+delta
443             call chainbuild_cart
444             uzap2=0.0d0
445             do ii=1,nfrag
446              qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
447      &         ,idummy,idummy)
448                uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
449      &          qinfrag(ii,iset))
450             enddo
451             do ii=1,npair
452                kstart=ifrag(1,ipair(1,ii,iset),iset)
453                kend=ifrag(2,ipair(1,ii,iset),iset)
454                lstart=ifrag(1,ipair(2,ii,iset),iset)
455                lend=ifrag(2,ipair(2,ii,iset),iset)
456                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
457                uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),
458      &          qinpair(ii,iset))
459             enddo
460             dc(j,i+nres)=cdummy(j,i)
461             call chainbuild_cart
462             uzap1=0.0d0
463              do ii=1,nfrag
464                qfrag(ii)=qwolynes(ifrag(1,ii,iset),
465      &          ifrag(2,ii,iset),.true.,idummy,idummy)
466                uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),
467      &          qinfrag(ii,iset))
468             enddo
469             do ii=1,npair
470                kstart=ifrag(1,ipair(1,ii,iset),iset)
471                kend=ifrag(2,ipair(1,ii,iset),iset)
472                lstart=ifrag(1,ipair(2,ii,iset),iset)
473                lend=ifrag(2,ipair(2,ii,iset),iset)
474                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
475                uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),
476      &          qinpair(ii,iset))
477             enddo
478             duxcartan(j,i)=(uzap2-uzap1)/(delta)            
479          enddo
480       enddo    
481       write(iout,*) "Numerical dUconst/ddc backbone "
482       do ii=0,nres
483         write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
484       enddo
485 c      write(iout,*) "Numerical dUconst/ddx side-chain "
486 c      do ii=1,nres
487 c         write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
488 c      enddo 
489       return
490       end
491 c---------------------------------------------------------------------------