changes to avoid PMPI_Scatterv: Invalid buffer pointer in Intel MPI
[unres.git] / source / unres / src_MD-M / 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.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),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 c-----------------------------------------------------------------------
253       subroutine dEconstrQ_num
254 c Calculating numerical dUconst/ddc and dUconst/ddx      
255       implicit real*8 (a-h,o-z)
256       include 'DIMENSIONS'
257       include 'COMMON.CONTROL'
258       include 'COMMON.VAR'
259       include 'COMMON.MD'
260 #ifndef LANG0
261       include 'COMMON.LANGEVIN'
262 #else
263       include 'COMMON.LANGEVIN.lang0'
264 #endif
265       include 'COMMON.CHAIN'
266       include 'COMMON.DERIV'
267       include 'COMMON.GEO'
268       include 'COMMON.LOCAL'
269       include 'COMMON.INTERACT'
270       include 'COMMON.IOUNITS'
271       include 'COMMON.NAMES'
272       include 'COMMON.TIME1'
273       double precision uzap1,uzap2
274       double precision dUcartan(3,0:MAXRES)
275      & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES)
276       integer kstart,kend,lstart,lend,idummy
277       double precision delta /1.0d-7/
278 c     For the backbone
279       do i=0,nres-1
280          do j=1,3
281             dUcartan(j,i)=0.0d0
282             cdummy(j,i)=dc(j,i)
283             dc(j,i)=dc(j,i)+delta
284             call chainbuild_cart
285             uzap2=0.0d0
286             do ii=1,nfrag
287                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
288      &         .true.,idummy,idummy)
289                uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
290      &            qinfrag(ii,iset))
291             enddo
292             do ii=1,npair
293                kstart=ifrag(1,ipair(1,ii,iset),iset)
294                kend=ifrag(2,ipair(1,ii,iset),iset)
295                lstart=ifrag(1,ipair(2,ii,iset),iset)
296                lend=ifrag(2,ipair(2,ii,iset),iset)
297                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
298                uzap2=uzap2+wpair(ii,iset)*
299      &             harmonic(qpair(ii),qinpair(ii,iset))
300             enddo
301             dc(j,i)=cdummy(j,i)
302             call chainbuild_cart
303             uzap1=0.0d0
304              do ii=1,nfrag
305                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
306      &         .true.,idummy,idummy)
307                uzap1=uzap1+wfrag(ii,iset)*
308      &                 harmonic(qfrag(ii),qinfrag(ii,iset))
309             enddo
310             do ii=1,npair
311                kstart=ifrag(1,ipair(1,ii,iset),iset)
312                kend=ifrag(2,ipair(1,ii,iset),iset)
313                lstart=ifrag(1,ipair(2,ii,iset),iset)
314                lend=ifrag(2,ipair(2,ii,iset),iset)
315                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
316                uzap1=uzap1+wpair(ii,iset)*
317      &            harmonic(qpair(ii),qinpair(ii,iset))
318             enddo
319             ducartan(j,i)=(uzap2-uzap1)/(delta)
320          enddo
321       enddo
322 c Calculating numerical gradients for dU/ddx
323       do i=0,nres-1
324          do j=1,3
325            duxcartan(j,i)=0.0d0
326          enddo
327          do j=1,3
328             cdummy(j,i)=dc(j,i+nres)
329             dc(j,i+nres)=dc(j,i+nres)+delta
330             call chainbuild_cart
331             uzap2=0.0d0
332             do ii=1,nfrag
333                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
334      &         .true.,idummy,idummy)
335                uzap2=uzap2+wfrag(ii,iset)*
336      &            harmonic(qfrag(ii),qinfrag(ii,iset))
337             enddo
338             do ii=1,npair
339                kstart=ifrag(1,ipair(1,ii,iset),iset)
340                kend=ifrag(2,ipair(1,ii,iset),iset)
341                lstart=ifrag(1,ipair(2,ii,iset),iset)
342                lend=ifrag(2,ipair(2,ii,iset),iset)
343                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
344                uzap2=uzap2+wpair(ii,iset)*
345      &             harmonic(qpair(ii),qinpair(ii,iset))
346             enddo
347             dc(j,i+nres)=cdummy(j,i)
348             call chainbuild_cart
349             uzap1=0.0d0
350              do ii=1,nfrag
351                qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
352      &         .true.,idummy,idummy)
353                uzap1=uzap1+wfrag(ii,iset)*
354      &            harmonic(qfrag(ii),qinfrag(ii,iset))
355             enddo
356             do ii=1,npair
357                kstart=ifrag(1,ipair(1,ii,iset),iset)
358                kend=ifrag(2,ipair(1,ii,iset),iset)
359                lstart=ifrag(1,ipair(2,ii,iset),iset)
360                lend=ifrag(2,ipair(2,ii,iset),iset)
361                qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
362                uzap1=uzap1+wpair(ii,iset)*
363      &               harmonic(qpair(ii),qinpair(ii,iset))
364             enddo
365             duxcartan(j,i)=(uzap2-uzap1)/(delta)
366          enddo
367       enddo    
368       write(iout,*) "Numerical dUconst/ddc backbone "
369       do ii=0,nres
370         write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
371       enddo
372       write(iout,*) "Numerical dUconst/ddx side-chain "
373       do ii=1,nres
374          write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
375       enddo 
376       return
377       end
378 c--------------------------------------------------------------------------- 
379       double precision function qcontrib(il,jl,il1,jl1)
380       implicit none
381       include 'DIMENSIONS'
382       include 'COMMON.IOUNITS'
383       include 'COMMON.CHAIN'
384       include 'COMMON.INTERACT'
385       include 'COMMON.MD'
386       integer i,j,k,il,jl,il1,jl1,nd
387       double precision dist
388       external dist
389       double precision dij1,dij2,dij3,dij4,d0ij1,d0ij2,d0ij3,d0ij4,fac,
390      &  fac1,ddave,ssij,ddqij
391       logical lprn /.false./
392       d0ij1=dsqrt((cref(1,jl)-cref(1,il))**2+
393      &           (cref(2,jl)-cref(2,il))**2+
394      &           (cref(3,jl)-cref(3,il))**2)
395       dij1=dist(il,jl)
396       ddave=(dij1-d0ij1)**2
397       nd=1
398       if (jl1.ne.jl) then
399         d0ij2=dsqrt((cref(1,jl1)-cref(1,il))**2+
400      &           (cref(2,jl1)-cref(2,il))**2+
401      &           (cref(3,jl1)-cref(3,il))**2)
402         dij2=dist(il,jl1)
403         ddave=ddave+(dij2-d0ij2)**2
404         nd=nd+1
405       endif
406       if (il1.ne.il) then
407         d0ij3=dsqrt((cref(1,jl)-cref(1,il1))**2+
408      &           (cref(2,jl)-cref(2,il1))**2+
409      &           (cref(3,jl)-cref(3,il1))**2)
410         dij3=dist(il1,jl)
411         ddave=ddave+(dij3-d0ij3)**2
412         nd=nd+1
413       endif
414       if (il1.ne.il .and. jl1.ne.jl) then
415         d0ij4=dsqrt((cref(1,jl1)-cref(1,il1))**2+
416      &           (cref(2,jl1)-cref(2,il1))**2+
417      &           (cref(3,jl1)-cref(3,il1))**2)
418         dij4=dist(il1,jl1)
419         ddave=ddave+(dij4-d0ij4)**2
420         nd=nd+1
421       endif
422       ddave=ddave/nd
423       if (lprn) then
424         write (iout,*) "il",il," jl",jl,
425      &  " itype",itype(il),itype(jl)," nd",nd
426         write (iout,*)"d0ij",d0ij1,d0ij2,d0ij3,d0ij4,
427      &  " dij",dij1,dij2,dij3,dij4," ddave",ddave
428         call flush(iout)
429       endif
430 c      ssij = (0.25d0*d0ij1)**2
431       if (il.ne.il1 .and. jl.ne.jl1) then
432         ssij = 16.0d0/(d0ij1*d0ij4)
433       else
434         ssij = 16.0d0/(d0ij1*d0ij1)
435       endif
436       qcontrib = dexp(-0.5d0*ddave*ssij)
437 c Compute gradient
438       fac1 = qcontrib*ssij/nd
439       fac = fac1*(dij1-d0ij1)/dij1
440       do k=1,3
441         ddqij = (c(k,il)-c(k,jl))*fac
442         dqwol(k,il)=dqwol(k,il)+ddqij
443         dqwol(k,jl)=dqwol(k,jl)-ddqij
444       enddo
445       if (jl1.ne.jl) then
446         fac = fac1*(dij2-d0ij2)/dij2
447         do k=1,3
448           ddqij = (c(k,il)-c(k,jl1))*fac
449           dqwol(k,il)=dqwol(k,il)+ddqij
450           dxqwol(k,jl)=dxqwol(k,jl)-ddqij
451         enddo
452       endif
453       if (il1.ne.il) then
454         fac = fac1*(dij3-d0ij3)/dij3
455         do k=1,3
456           ddqij = (c(k,il1)-c(k,jl))*fac
457           dxqwol(k,il)=dxqwol(k,il)+ddqij
458           dqwol(k,jl)=dqwol(k,jl)-ddqij
459         enddo
460       endif
461       if (il1.ne.il .and. jl1.ne.jl) then
462         fac = fac1*(dij4-d0ij4)/dij4
463         do k=1,3
464           ddqij = (c(k,il1)-c(k,jl1))*fac
465           dxqwol(k,il)=dxqwol(k,il)+ddqij
466           dxqwol(k,jl)=dxqwol(k,jl)-ddqij
467         enddo
468       endif
469       return
470       end