changes to avoid PMPI_Scatterv: Invalid buffer pointer in Intel MPI
[unres.git] / source / unres / src_MD-M / q_measure3.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.MD'
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,iset),ifrag(2,i,iset),.true.,idummy,idummy)
199       enddo
200 c      stop
201       do i=1,npair
202          kstart=ifrag(1,ipair(1,i,iset),iset)
203          kend=ifrag(2,ipair(1,i,iset),iset)
204          lstart=ifrag(1,ipair(2,i,iset),iset)
205          lend=ifrag(2,ipair(2,i,iset),iset)
206          qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
207          Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
208 c  Calculating dU/dQ
209          Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
210 c Calculating dQ/dXi
211          do ii=0,nres
212             do j=1,3
213                duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
214                dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
215             enddo
216          enddo
217       enddo
218 c      write(iout,*) "Uconst inside subroutine ", Uconst
219 c Transforming the gradients from Cs to dCs for the backbone
220       do i=0,nres
221          do j=i+1,nres
222            do k=1,3
223              dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
224            enddo
225          enddo
226       enddo
227 c  Transforming the gradients from Cs to dCs for the side chains      
228       do i=1,nres
229          do j=1,3
230            dudxconst(j,i)=duxconst(j,i)
231          enddo
232       enddo
233 c      write(iout,*) "dU/dc backbone "
234 c       do ii=0,nres
235 c        write(iout,'(i5,3e15.5)') ii, (duconst(j,ii),j=1,3)
236 c      enddo      
237 c      write(iout,*) "dU/dX side chain "
238 c      do ii=1,nres
239 c            write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
240 c      enddo
241 c      write(iout,*) "dU/ddc backbone "
242 c       do ii=0,nres
243 c        write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
244 c      enddo      
245 c      write(iout,*) "dU/ddX side chain "
246 c      do ii=1,nres
247 c            write(iout,'(i5,3e15.5)') ii,(dudxconst(j,ii),j=1,3)
248 c      enddo
249 c Calculating numerical gradients of dUconst/ddc and dUconst/ddx
250 c      call dEconstrQ_num      
251       return
252       end
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)*
291      &                harmonic(qfrag(ii),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 c--------------------------------------------------------------------------- 
380       double precision function qcontrib(il,jl,il1,jl1)
381       implicit none
382       include 'DIMENSIONS'
383       include 'COMMON.IOUNITS'
384       include 'COMMON.CHAIN'
385       include 'COMMON.INTERACT'
386       include 'COMMON.MD'
387       include 'COMMON.LOCAL'
388       integer i,j,k,il,jl,il1,jl1,nd,itl,jtl
389       double precision dist
390       external dist
391       double precision dij,dij1,d0ij,d0ij1,om1,om2,om12,om10,om20,om120
392      &  ,fac,fac1,ddave,ssij,ddqij,d0ii1,d0jj1,rij,eom1,eom2,eom12
393       double precision u(3),v(3),er(3),er0(3),dcosom1(3),dcosom2(3),
394      &  aux1,aux2
395       double precision scalar
396       external scalar
397       logical lprn /.false./
398       if (lprn) write (iout,*) "il",il," jl",jl," il1",il1," jl1",jl1
399       d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
400      &           (cref(2,jl)-cref(2,il))**2+
401      &           (cref(3,jl)-cref(3,il))**2)
402       dij=dist(il,jl)
403       dij1=dist(il1,jl1)
404       do i=1,3
405         er(i)=(c(i,jl1)-c(i,il1))/dij1
406       enddo
407       do i=1,3
408         er0(i)=cref(i,jl1)-cref(i,il1)
409       enddo
410       d0ij1=dsqrt(scalar(er0,er0))
411       do i=1,3
412         er0(i)=er0(i)/d0ij1
413       enddo
414       if (il.ne.il1 .or. jl.ne.jl1) then
415         ddave=0.5d0*((dij-d0ij)**2+(dij1-d0ij1)**2)
416         nd=2
417       else
418         ddave=(dij-d0ij)**2
419         nd=1
420       endif
421       if (il.ne.il1) then
422         do i=1,3
423           u(i)=cref(i,il1)-cref(i,il)
424         enddo
425         d0ii1=dsqrt(scalar(u,u))
426         do i=1,3
427           u(i)=u(i)/d0ii1
428         enddo
429         if (lprn) then
430         write (iout,*) "u",(u(i),i=1,3)
431         write (iout,*) "er0",(er0(i),i=1,3)
432         om10=scalar(er0,u)
433         om1=scalar(er,dc_norm(1,il1))
434         write (iout,*) "om10",om10," om1",om1
435         endif
436       else
437         om1=0.0d0
438         om10=0.0d0
439       endif
440       if (jl.ne.jl1) then
441         do i=1,3
442           v(i)=cref(i,jl1)-cref(i,jl)
443         enddo
444         d0jj1=dsqrt(scalar(v,v))
445         do i=1,3
446           v(i)=v(i)/d0jj1
447         enddo
448         if (lprn) then
449         write (iout,*) "v",(v(i),i=1,3)
450         write (iout,*) "er0",(er0(i),i=1,3)
451         om20=scalar(er,v)
452         om2=scalar(er,dc_norm(1,jl1))
453         write (iout,*) "om20",om20," om2",om2
454         endif
455       else
456         om2=0.0d0
457         om20=0.0d0
458       endif
459       if (il.ne.il1 .and. jl.ne.jl1) then
460         om120=scalar(u,v)
461         om12=scalar(dc_norm(1,il1),dc_norm(1,jl1))
462       else
463         om12=0.0d0
464         om120=0.0d0
465       endif
466       if (lprn) then
467         write (iout,*) "il",il," jl",jl,itype(il),itype(jl)
468         write (iout,*)"d0ij",d0ij," om10",om10," om20",om20,
469      &   " om120",om120,
470      &  " dij",dij," om1",om1," om2",om2," om12",om12
471         call flush(iout)
472       endif
473       ssij = 16.0d0/(d0ij*d0ij)
474       qcontrib = dexp(-0.5d0*(ddave*ssij+((om1-om10)**2
475      &                       +(om2-om20)**2+(om12-om120)**2)))
476       if (lprn) write (iout,*) "ssij",ssij," qcontrib",qcontrib
477 c      qcontrib = dexp(-0.5d0*(ddave*ssij)+(om1-om10)**2+(om2-om20)**2)
478 c      qcontrib = dexp(-0.5d0*(ddave*ssij))
479 c Compute gradient - radial component
480       fac1 = qcontrib*ssij/nd
481       fac = fac1*(dij-d0ij)/dij
482       do k=1,3
483         ddqij = (c(k,il)-c(k,jl))*fac
484         dqwol(k,il)=dqwol(k,il)+ddqij
485         dqwol(k,jl)=dqwol(k,jl)-ddqij
486       enddo
487       if (il1.ne.il .or. jl1.ne.jl) then
488         fac = fac1*(dij1-d0ij1)/dij1
489         do k=1,3
490           ddqij = (c(k,il1)-c(k,jl1))*fac
491           if (il1.ne.il) then
492             dxqwol(k,il)=dxqwol(k,il)+ddqij
493           else
494             dqwol(k,il)=dqwol(k,il)+ddqij
495           endif
496           if (jl1.ne.jl) then
497             dxqwol(k,jl)=dxqwol(k,jl)-ddqij
498           else
499             dqwol(k,jl)=dqwol(k,jl)-ddqij
500           endif
501         enddo
502       endif
503 c      return
504 c Orientational contributions
505       rij=1.0d0/dij1
506       eom1=qcontrib*(om1-om10)
507       eom2=qcontrib*(om2-om20)
508       eom12=qcontrib*(om12-om120)
509       do k=1,3
510         dcosom1(k)=rij*(dc_norm(k,il1)-om1*er(k))
511         dcosom2(k)=rij*(dc_norm(k,jl1)-om2*er(k))
512       enddo
513       do k=1,3
514         ddqij=eom1*dcosom1(k)+eom2*dcosom2(k)
515         aux1=(eom12*(dc_norm(k,jl1)-om12*dc_norm(k,il1))
516      &            +eom1*(er(k)-om1*dc_norm(k,il1)))*vbld_inv(il1)
517         aux2=(eom12*(dc_norm(k,il1)-om12*dc_norm(k,jl1))
518      &            +eom2*(er(k)-om2*dc_norm(k,jl1)))*vbld_inv(jl1)
519         dqwol(k,il)=dqwol(k,il)-ddqij-aux1
520         dqwol(k,jl)=dqwol(k,jl)+ddqij-aux2
521         dxqwol(k,il)=dxqwol(k,il)-ddqij+aux1
522 c     &            +(eom12*(dc_norm(k,jl1)-om12*dc_norm(k,il1))
523 c     &            +eom1*(er(k)-om1*dc_norm(k,il1)))*vbld_inv(il1)
524         dxqwol(k,jl)=dxqwol(k,jl)+ddqij+aux2
525 c     &            +(eom12*(dc_norm(k,il1)-om12*dc_norm(k,jl1))
526 c     &            +eom2*(er(k)-om2*dc_norm(k,jl1)))*vbld_inv(jl1)
527       enddo
528       return
529       end