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