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