energy_dec for esccor and some other debug printout
[unres.git] / source / unres / src_CSA / checkder_p.F
1       subroutine check_cartgrad
2 C Check the gradient of Cartesian coordinates in internal coordinates.
3       implicit real*8 (a-h,o-z)
4       include 'DIMENSIONS'
5       include 'COMMON.IOUNITS'
6       include 'COMMON.VAR'
7       include 'COMMON.CHAIN'
8       include 'COMMON.GEO'
9       include 'COMMON.LOCAL'
10       include 'COMMON.DERIV'
11       dimension temp(6,maxres),xx(3),gg(3)
12       indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
13 *
14 * Check the gradient of the virtual-bond and SC vectors in the internal
15 * coordinates.
16 *    
17       aincr=1.0d-7  
18       aincr2=5.0d-8   
19       call cartder
20       write (iout,'(a)') '**************** dx/dalpha'
21       write (iout,'(a)')
22       do i=2,nres-1
23         alphi=alph(i)
24         alph(i)=alph(i)+aincr
25         do k=1,3
26           temp(k,i)=dc(k,nres+i)
27         enddo
28         call chainbuild
29         do k=1,3
30           gg(k)=(dc(k,nres+i)-temp(k,i))/aincr
31           xx(k)=dabs((gg(k)-dxds(k,i))/(aincr*dabs(dxds(k,i))+aincr))
32         enddo
33         write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)') 
34      &  i,(gg(k),k=1,3),(dxds(k,i),k=1,3),(xx(k),k=1,3)
35         write (iout,'(a)')
36         alph(i)=alphi
37         call chainbuild
38       enddo
39       write (iout,'(a)')
40       write (iout,'(a)') '**************** dx/domega'
41       write (iout,'(a)')
42       do i=2,nres-1
43         omegi=omeg(i)
44         omeg(i)=omeg(i)+aincr
45         do k=1,3
46           temp(k,i)=dc(k,nres+i)
47         enddo
48         call chainbuild
49         do k=1,3
50           gg(k)=(dc(k,nres+i)-temp(k,i))/aincr
51           xx(k)=dabs((gg(k)-dxds(k+3,i))/
52      &          (aincr*dabs(dxds(k+3,i))+aincr))
53         enddo
54         write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)') 
55      &      i,(gg(k),k=1,3),(dxds(k+3,i),k=1,3),(xx(k),k=1,3)
56         write (iout,'(a)')
57         omeg(i)=omegi
58         call chainbuild
59       enddo
60       write (iout,'(a)')
61       write (iout,'(a)') '**************** dx/dtheta'
62       write (iout,'(a)')
63       do i=3,nres
64         theti=theta(i)
65         theta(i)=theta(i)+aincr
66         do j=i-1,nres-1
67           do k=1,3
68             temp(k,j)=dc(k,nres+j)
69           enddo
70         enddo
71         call chainbuild
72         do j=i-1,nres-1
73           ii = indmat(i-2,j)
74 c         print *,'i=',i-2,' j=',j-1,' ii=',ii
75           do k=1,3
76             gg(k)=(dc(k,nres+j)-temp(k,j))/aincr
77             xx(k)=dabs((gg(k)-dxdv(k,ii))/
78      &            (aincr*dabs(dxdv(k,ii))+aincr))
79           enddo
80           write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') 
81      &        i,j,(gg(k),k=1,3),(dxdv(k,ii),k=1,3),(xx(k),k=1,3)
82           write(iout,'(a)')
83         enddo
84         write (iout,'(a)')
85         theta(i)=theti
86         call chainbuild
87       enddo
88       write (iout,'(a)') '***************** dx/dphi'
89       write (iout,'(a)')
90       do i=4,nres
91         phi(i)=phi(i)+aincr
92         do j=i-1,nres-1
93           do k=1,3
94             temp(k,j)=dc(k,nres+j)
95           enddo
96         enddo
97         call chainbuild
98         do j=i-1,nres-1
99           ii = indmat(i-2,j)
100 c         print *,'ii=',ii
101           do k=1,3
102             gg(k)=(dc(k,nres+j)-temp(k,j))/aincr
103             xx(k)=dabs((gg(k)-dxdv(k+3,ii))/
104      &            (aincr*dabs(dxdv(k+3,ii))+aincr))
105           enddo
106           write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') 
107      &        i,j,(gg(k),k=1,3),(dxdv(k+3,ii),k=1,3),(xx(k),k=1,3)
108           write(iout,'(a)')
109         enddo
110         phi(i)=phi(i)-aincr
111         call chainbuild
112       enddo
113       write (iout,'(a)') '****************** ddc/dtheta'
114       do i=1,nres-2
115         thet=theta(i+2)
116         theta(i+2)=thet+aincr
117         do j=i,nres
118           do k=1,3 
119             temp(k,j)=dc(k,j)
120           enddo
121         enddo
122         call chainbuild 
123         do j=i+1,nres-1
124           ii = indmat(i,j)
125 c         print *,'ii=',ii
126           do k=1,3
127             gg(k)=(dc(k,j)-temp(k,j))/aincr
128             xx(k)=dabs((gg(k)-dcdv(k,ii))/
129      &           (aincr*dabs(dcdv(k,ii))+aincr))
130           enddo
131           write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') 
132      &           i,j,(gg(k),k=1,3),(dcdv(k,ii),k=1,3),(xx(k),k=1,3)
133           write (iout,'(a)')
134         enddo
135         do j=1,nres
136           do k=1,3
137             dc(k,j)=temp(k,j)
138           enddo 
139         enddo
140         theta(i+2)=thet
141       enddo    
142       write (iout,'(a)') '******************* ddc/dphi'
143       do i=1,nres-3
144         phii=phi(i+3)
145         phi(i+3)=phii+aincr
146         do j=1,nres
147           do k=1,3 
148             temp(k,j)=dc(k,j)
149           enddo
150         enddo
151         call chainbuild 
152         do j=i+2,nres-1
153           ii = indmat(i+1,j)
154 c         print *,'ii=',ii
155           do k=1,3
156             gg(k)=(dc(k,j)-temp(k,j))/aincr
157             xx(k)=dabs((gg(k)-dcdv(k+3,ii))/
158      &           (aincr*dabs(dcdv(k+3,ii))+aincr))
159           enddo
160           write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') 
161      &         i,j,(gg(k),k=1,3),(dcdv(k+3,ii),k=1,3),(xx(k),k=1,3)
162           write (iout,'(a)')
163         enddo
164         do j=1,nres
165           do k=1,3
166             dc(k,j)=temp(k,j)
167           enddo
168         enddo
169         phi(i+3)=phii   
170       enddo   
171       return
172       end
173 C----------------------------------------------------------------------------
174       subroutine check_ecart
175 C Check the gradient of the energy in Cartesian coordinates. 
176       implicit real*8 (a-h,o-z)
177       include 'DIMENSIONS'
178       include 'COMMON.CHAIN'
179       include 'COMMON.DERIV'
180       include 'COMMON.IOUNITS'
181       include 'COMMON.VAR'
182       include 'COMMON.CONTACTS'
183 #ifdef MOMENT
184       include 'COMMON.CONTACTS.MOMENT'
185 #endif  
186       common /srutu/ icall
187       dimension ggg(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),g(maxvar)
188       dimension grad_s(6,maxres)
189       double precision energia(0:n_ene),energia1(0:n_ene)
190       integer uiparm(1)
191       double precision urparm(1)
192       external fdum
193       icg=1
194       nf=0
195       nfl=0                
196       call zerograd
197       aincr=1.0D-7
198       print '(a)','CG processor',me,' calling CHECK_CART.'
199       nf=0
200       icall=0
201       call geom_to_var(nvar,x)
202       call etotal(energia(0))
203       etot=energia(0)
204       call enerprint(energia(0))
205       call gradient(nvar,x,nf,g,uiparm,urparm,fdum)
206       icall =1
207       do i=1,nres
208         write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
209       enddo
210       do i=1,nres
211         do j=1,3
212           grad_s(j,i)=gradc(j,i,icg)
213           grad_s(j+3,i)=gradx(j,i,icg)
214         enddo
215       enddo
216       call flush(iout)
217       write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
218       do i=1,nres
219         do j=1,3
220           xx(j)=c(j,i+nres)
221           ddc(j)=dc(j,i) 
222           ddx(j)=dc(j,i+nres)
223         enddo
224         do j=1,3
225           dc(j,i)=dc(j,i)+aincr
226           do k=i+1,nres
227             c(j,k)=c(j,k)+aincr
228             c(j,k+nres)=c(j,k+nres)+aincr
229           enddo
230           call etotal(energia1(0))
231           etot1=energia1(0)
232           ggg(j)=(etot1-etot)/aincr
233           dc(j,i)=ddc(j)
234           do k=i+1,nres
235             c(j,k)=c(j,k)-aincr
236             c(j,k+nres)=c(j,k+nres)-aincr
237           enddo
238         enddo
239         do j=1,3
240           c(j,i+nres)=c(j,i+nres)+aincr
241           dc(j,i+nres)=dc(j,i+nres)+aincr
242           call etotal(energia1(0))
243           etot1=energia1(0)
244           ggg(j+3)=(etot1-etot)/aincr
245           c(j,i+nres)=xx(j)
246           dc(j,i+nres)=ddx(j)
247         enddo
248         write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/)')
249      &   i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6)
250       enddo
251       return
252       end
253 c----------------------------------------------------------------------------
254       subroutine check_ecartint
255 C Check the gradient of the energy in Cartesian coordinates. 
256       implicit real*8 (a-h,o-z)
257       include 'DIMENSIONS'
258       include 'COMMON.CONTROL'
259       include 'COMMON.CHAIN'
260       include 'COMMON.DERIV'
261       include 'COMMON.IOUNITS'
262       include 'COMMON.VAR'
263       include 'COMMON.CONTACTS'
264 #ifdef MOMENT
265       include 'COMMON.CONTACTS.MOMENT'
266 #endif  
267       include 'COMMON.MD_'
268       include 'COMMON.LOCAL'
269       include 'COMMON.SPLITELE'
270       common /srutu/ icall
271       dimension ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),
272      &  g(maxvar)
273       dimension dcnorm_safe(3),dxnorm_safe(3)
274       dimension grad_s(6,0:maxres),grad_s1(6,0:maxres)
275       double precision phi_temp(maxres),theta_temp(maxres),
276      &  alph_temp(maxres),omeg_temp(maxres)
277       double precision energia(0:n_ene),energia1(0:n_ene)
278       integer uiparm(1)
279       double precision urparm(1)
280       external fdum
281       r_cut=2.0d0
282       rlambd=0.3d0
283       icg=1
284       nf=0
285       nfl=0                
286       call intout
287 c      call intcartderiv
288 c      call checkintcartgrad
289       call zerograd
290       aincr=1.0D-5
291       write(iout,*) 'Calling CHECK_ECARTINT.'
292       nf=0
293       icall=0
294       call geom_to_var(nvar,x)
295       if (.not.split_ene) then
296         call etotal(energia(0))
297         etot=energia(0)
298         call enerprint(energia(0))
299         call flush(iout)
300         write (iout,*) "enter cartgrad"
301         call flush(iout)
302         call cartgrad
303         write (iout,*) "exit cartgrad"
304         call flush(iout)
305         icall =1
306         do i=1,nres
307           write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
308         enddo
309         do j=1,3
310           grad_s(j,0)=gcart(j,0)
311         enddo
312         do i=1,nres
313           do j=1,3
314             grad_s(j,i)=gcart(j,i)
315             grad_s(j+3,i)=gxcart(j,i)
316           enddo
317         enddo
318       else
319 !- split gradient check
320          write (iout,*) "split_ene not supported"
321 c        call zerograd
322 c        call etotal_long(energia(0))
323 c        call enerprint(energia(0))
324 c        call flush(iout)
325 c        write (iout,*) "enter cartgrad"
326 c        call flush(iout)
327 c        call cartgrad
328 c        write (iout,*) "exit cartgrad"
329 c        call flush(iout)
330 c        icall =1
331 c        write (iout,*) "longrange grad"
332 c        do i=1,nres
333 c          write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
334 c     &    (gxcart(j,i),j=1,3)
335 c        enddo
336 c        do j=1,3
337 c          grad_s(j,0)=gcart(j,0)
338 c        enddo
339 c        do i=1,nres
340 c          do j=1,3
341 c            grad_s(j,i)=gcart(j,i)
342 c            grad_s(j+3,i)=gxcart(j,i)
343 c          enddo
344 c        enddo
345 c        call zerograd
346 c        call etotal_short(energia(0))
347 c        call enerprint(energia(0))
348 c        call flush(iout)
349 c        write (iout,*) "enter cartgrad"
350 c        call flush(iout)
351 c        call cartgrad
352 c        write (iout,*) "exit cartgrad"
353 c        call flush(iout)
354 c        icall =1
355 c        write (iout,*) "shortrange grad"
356 c        do i=1,nres
357 c          write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
358 c     &    (gxcart(j,i),j=1,3)
359 c        enddo
360 c        do j=1,3
361 c          grad_s1(j,0)=gcart(j,0)
362 c        enddo
363 c        do i=1,nres
364 c          do j=1,3
365 c            grad_s1(j,i)=gcart(j,i)
366 c            grad_s1(j+3,i)=gxcart(j,i)
367 c          enddo
368 c        enddo
369       endif
370       write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
371       do i=0,nres
372         do j=1,3
373           xx(j)=c(j,i+nres)
374           ddc(j)=dc(j,i) 
375           ddx(j)=dc(j,i+nres)
376           do k=1,3
377             dcnorm_safe(k)=dc_norm(k,i)
378             dxnorm_safe(k)=dc_norm(k,i+nres)
379           enddo
380         enddo
381         do j=1,3
382           dc(j,i)=ddc(j)+aincr
383           call chainbuild_cart
384 #ifdef MPI
385 c Broadcast the order to compute internal coordinates to the slaves.
386 c          if (nfgtasks.gt.1)
387 c     &      call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
388 #endif
389 c          call int_from_cart1(.false.)
390           if (.not.split_ene) then
391             call etotal(energia1(0))
392             etot1=energia1(0)
393           else
394 !- split gradient
395 c            call etotal_long(energia1(0))
396 c            etot11=energia1(0)
397 c            call etotal_short(energia1(0))
398 c            etot12=energia1(0)
399 c            write (iout,*) "etot11",etot11," etot12",etot12
400           endif
401 !- end split gradient
402 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
403           dc(j,i)=ddc(j)-aincr
404           call chainbuild_cart
405 c          call int_from_cart1(.false.)
406           if (.not.split_ene) then
407             call etotal(energia1(0))
408             etot2=energia1(0)
409             ggg(j)=(etot1-etot2)/(2*aincr)
410           else
411 !- split gradient
412 c            call etotal_long(energia1(0))
413 c            etot21=energia1(0)
414 c           ggg(j)=(etot11-etot21)/(2*aincr)
415 c            call etotal_short(energia1(0))
416 c            etot22=energia1(0)
417 c           ggg1(j)=(etot12-etot22)/(2*aincr)
418 !- end split gradient
419 c            write (iout,*) "etot21",etot21," etot22",etot22
420           endif
421 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
422           dc(j,i)=ddc(j)
423           call chainbuild_cart
424         enddo
425         do j=1,3
426           dc(j,i+nres)=ddx(j)+aincr
427           call chainbuild_cart
428 c          write (iout,*) "i",i," j",j," dxnorm+ and dxnorm"
429 c          write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
430 c          write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
431 c          write (iout,*) "dxnormnorm",dsqrt(
432 c     &  dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
433 c          write (iout,*) "dxnormnormsafe",dsqrt(
434 c     &      dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
435 c          write (iout,*)
436           if (.not.split_ene) then
437             call etotal(energia1(0))
438             etot1=energia1(0)
439           else
440 !- split gradient
441 c            call etotal_long(energia1(0))
442 c            etot11=energia1(0)
443 c            call etotal_short(energia1(0))
444 c            etot12=energia1(0)
445           endif
446 !- end split gradient
447 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
448           dc(j,i+nres)=ddx(j)-aincr
449           call chainbuild_cart
450 c          write (iout,*) "i",i," j",j," dxnorm- and dxnorm"
451 c          write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
452 c          write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
453 c          write (iout,*) 
454 c          write (iout,*) "dxnormnorm",dsqrt(
455 c     &  dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
456 c          write (iout,*) "dxnormnormsafe",dsqrt(
457 c     &      dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
458           if (.not.split_ene) then
459             call etotal(energia1(0))
460             etot2=energia1(0)
461             ggg(j+3)=(etot1-etot2)/(2*aincr)
462           else
463 !- split gradient
464 c            call etotal_long(energia1(0))
465 c            etot21=energia1(0)
466 c           ggg(j+3)=(etot11-etot21)/(2*aincr)
467 c           call etotal_short(energia1(0))
468 c            etot22=energia1(0)
469 c           ggg1(j+3)=(etot12-etot22)/(2*aincr)
470 !- end split gradient
471           endif
472 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
473           dc(j,i+nres)=ddx(j)
474           call chainbuild_cart
475         enddo
476         write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
477      &   i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
478         if (split_ene) then
479           write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
480      &   i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
481      &   k=1,6)
482          write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
483      &   i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
484      &   ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
485         endif
486       enddo
487       return
488       end
489 c-------------------------------------------------------------------------
490       subroutine int_from_cart1(lprn)
491       implicit real*8 (a-h,o-z)
492       include 'DIMENSIONS'
493 #ifdef MPI
494       include 'mpif.h'
495       integer ierror
496 #endif
497       include 'COMMON.IOUNITS'
498       include 'COMMON.VAR'
499       include 'COMMON.CHAIN'
500       include 'COMMON.GEO'
501       include 'COMMON.INTERACT'
502       include 'COMMON.LOCAL'
503       include 'COMMON.NAMES'
504       include 'COMMON.SETUP'
505       include 'COMMON.TIME1'
506       logical lprn 
507       if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
508 #ifdef TIMING
509       time01=MPI_Wtime()
510 #endif
511 #if defined(PARINT) && defined(MPI)
512       do i=iint_start,iint_end
513 #else
514       do i=2,nres
515 #endif
516         dnorm1=dist(i-1,i)
517         dnorm2=dist(i,i+1) 
518         do j=1,3
519           c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1
520      &     +(c(j,i+1)-c(j,i))/dnorm2)
521         enddo
522         be=0.0D0
523         if (i.gt.2) then
524         if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
525         if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
526          tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
527         endif
528         if (itype(i-1).ne.10) then
529          tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
530          omicron(1,i)=alpha(i-2,i-1,i-1+nres)
531          omicron(2,i)=alpha(i-1+nres,i-1,i)
532         endif
533         if (itype(i).ne.10) then
534          tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
535         endif
536         endif
537         omeg(i)=beta(nres+i,i,maxres2,i+1)
538         alph(i)=alpha(nres+i,i,maxres2)
539         theta(i+1)=alpha(i-1,i,i+1)
540         vbld(i)=dist(i-1,i)
541         vbld_inv(i)=1.0d0/vbld(i)
542         vbld(nres+i)=dist(nres+i,i)
543         if (itype(i).ne.10) then
544           vbld_inv(nres+i)=1.0d0/vbld(nres+i)
545         else
546           vbld_inv(nres+i)=0.0d0
547         endif
548       enddo   
549 #if defined(PARINT) && defined(MPI)
550        if (nfgtasks1.gt.1) then
551 cd       write(iout,*) "iint_start",iint_start," iint_count",
552 cd     &   (iint_count(i),i=0,nfgtasks-1)," iint_displ",
553 cd     &   (iint_displ(i),i=0,nfgtasks-1)
554 cd       write (iout,*) "Gather vbld backbone"
555 cd       call flush(iout)
556        time00=MPI_Wtime()
557        call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
558      &   MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
559      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
560 cd       write (iout,*) "Gather vbld_inv"
561 cd       call flush(iout)
562        call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
563      &   MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
564      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
565 cd       write (iout,*) "Gather vbld side chain"
566 cd       call flush(iout)
567        call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
568      &  MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
569      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
570 cd       write (iout,*) "Gather vbld_inv side chain"
571 cd       call flush(iout)
572        call MPI_Allgatherv(vbld_inv(iint_start+nres),
573      &   iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
574      &   iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
575 cd       write (iout,*) "Gather theta"
576 cd       call flush(iout)
577        call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
578      &   MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
579      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
580 cd       write (iout,*) "Gather phi"
581 cd       call flush(iout)
582        call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
583      &   MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
584      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
585 #ifdef CRYST_SC
586 cd       write (iout,*) "Gather alph"
587 cd       call flush(iout)
588        call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
589      &   MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
590      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
591 cd       write (iout,*) "Gather omeg"
592 cd       call flush(iout)
593        call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
594      &   MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
595      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
596 #endif
597        time_gather=time_gather+MPI_Wtime()-time00
598       endif
599 #endif
600       do i=1,nres-1
601         do j=1,3
602           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
603         enddo
604       enddo
605       do i=2,nres-1
606         do j=1,3
607           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
608         enddo
609       enddo
610       if (lprn) then
611       do i=2,nres
612        write (iout,1212) restyp(itype(i)),i,vbld(i),
613      &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
614      &rad2deg*alph(i),rad2deg*omeg(i)
615       enddo
616       endif
617  1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
618 #ifdef TIMING
619       time_intfcart=time_intfcart+MPI_Wtime()-time01
620 #endif
621       return
622       end
623 c----------------------------------------------------------------------------
624       subroutine check_eint
625 C Check the gradient of energy in internal coordinates.
626       implicit real*8 (a-h,o-z)
627       include 'DIMENSIONS'
628       include 'COMMON.CHAIN'
629       include 'COMMON.DERIV'
630       include 'COMMON.IOUNITS'
631       include 'COMMON.VAR'
632       include 'COMMON.GEO'
633       common /srutu/ icall
634       dimension x(maxvar),gana(maxvar),gg(maxvar)
635       integer uiparm(1)
636       double precision urparm(1)
637       double precision energia(0:n_ene),energia1(0:n_ene),
638      &  energia2(0:n_ene)
639       character*6 key
640       external fdum
641       call zerograd
642       aincr=1.0D-7
643       print '(a)','Calling CHECK_INT.'
644       nf=0
645       nfl=0
646       icg=1
647       call geom_to_var(nvar,x)
648       call var_to_geom(nvar,x)
649       call chainbuild
650       icall=1
651       print *,'ICG=',ICG
652       call etotal(energia(0))
653       etot = energia(0)
654       call enerprint(energia(0))
655       print *,'ICG=',ICG
656 #ifdef MPL
657       if (MyID.ne.BossID) then
658         call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
659         nf=x(nvar+1)
660         nfl=x(nvar+2)
661         icg=x(nvar+3)
662       endif
663 #endif
664       nf=1
665       nfl=3
666 cd    write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
667       call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
668 cd    write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
669       icall=1
670       do i=1,nvar
671         xi=x(i)
672         x(i)=xi-0.5D0*aincr
673         call var_to_geom(nvar,x)
674         call chainbuild
675         call etotal(energia1(0))
676         etot1=energia1(0)
677         x(i)=xi+0.5D0*aincr
678         call var_to_geom(nvar,x)
679         call chainbuild
680         call etotal(energia2(0))
681         etot2=energia2(0)
682         gg(i)=(etot2-etot1)/aincr
683         write (iout,*) i,etot1,etot2
684         x(i)=xi
685       enddo
686       write (iout,'(/2a)')' Variable        Numerical       Analytical',
687      &    '     RelDiff*100% '
688       do i=1,nvar
689         if (i.le.nphi) then
690           ii=i
691           key = ' phi'
692         else if (i.le.nphi+ntheta) then
693           ii=i-nphi
694           key=' theta'
695         else if (i.le.nphi+ntheta+nside) then
696            ii=i-(nphi+ntheta)
697            key=' alpha'
698         else 
699            ii=i-(nphi+ntheta+nside)
700            key=' omega'
701         endif
702         write (iout,'(i3,a,i3,3(1pd16.6))') 
703      & i,key,ii,gg(i),gana(i),
704      & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)
705       enddo
706       return
707       end