a50bc2a790756e0caaa5a15601a15055ed8dbff0
[unres.git] / 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       common /srutu/ icall
184       dimension ggg(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),g(maxvar)
185       dimension grad_s(6,maxres)
186       double precision energia(0:n_ene),energia1(0:n_ene)
187       integer uiparm(1)
188       double precision urparm(1)
189       external fdum
190       icg=1
191       nf=0
192       nfl=0                
193       call zerograd
194       aincr=1.0D-7
195       print '(a)','CG processor',me,' calling CHECK_CART.'
196       nf=0
197       icall=0
198       call geom_to_var(nvar,x)
199       call etotal(energia(0))
200       etot=energia(0)
201       call enerprint(energia(0))
202       call gradient(nvar,x,nf,g,uiparm,urparm,fdum)
203       icall =1
204       do i=1,nres
205         write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
206       enddo
207       do i=1,nres
208         do j=1,3
209           grad_s(j,i)=gradc(j,i,icg)
210           grad_s(j+3,i)=gradx(j,i,icg)
211         enddo
212       enddo
213       call flush(iout)
214       write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
215       do i=1,nres
216         do j=1,3
217           xx(j)=c(j,i+nres)
218           ddc(j)=dc(j,i) 
219           ddx(j)=dc(j,i+nres)
220         enddo
221         do j=1,3
222           dc(j,i)=dc(j,i)+aincr
223           do k=i+1,nres
224             c(j,k)=c(j,k)+aincr
225             c(j,k+nres)=c(j,k+nres)+aincr
226           enddo
227           call etotal(energia1(0))
228           etot1=energia1(0)
229           ggg(j)=(etot1-etot)/aincr
230           dc(j,i)=ddc(j)
231           do k=i+1,nres
232             c(j,k)=c(j,k)-aincr
233             c(j,k+nres)=c(j,k+nres)-aincr
234           enddo
235         enddo
236         do j=1,3
237           c(j,i+nres)=c(j,i+nres)+aincr
238           dc(j,i+nres)=dc(j,i+nres)+aincr
239           call etotal(energia1(0))
240           etot1=energia1(0)
241           ggg(j+3)=(etot1-etot)/aincr
242           c(j,i+nres)=xx(j)
243           dc(j,i+nres)=ddx(j)
244         enddo
245         write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/)')
246      &   i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6)
247       enddo
248       return
249       end
250 c----------------------------------------------------------------------------
251       subroutine check_ecartint
252 C Check the gradient of the energy in Cartesian coordinates. 
253       implicit real*8 (a-h,o-z)
254       include 'DIMENSIONS'
255       include 'COMMON.CONTROL'
256       include 'COMMON.CHAIN'
257       include 'COMMON.DERIV'
258       include 'COMMON.IOUNITS'
259       include 'COMMON.VAR'
260       include 'COMMON.CONTACTS'
261       include 'COMMON.MD_'
262       include 'COMMON.LOCAL'
263       include 'COMMON.SPLITELE'
264       common /srutu/ icall
265       dimension ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),
266      &  g(maxvar)
267       dimension dcnorm_safe(3),dxnorm_safe(3)
268       dimension grad_s(6,0:maxres),grad_s1(6,0:maxres)
269       double precision phi_temp(maxres),theta_temp(maxres),
270      &  alph_temp(maxres),omeg_temp(maxres)
271       double precision energia(0:n_ene),energia1(0:n_ene)
272       integer uiparm(1)
273       double precision urparm(1)
274       external fdum
275       r_cut=2.0d0
276       rlambd=0.3d0
277       icg=1
278       nf=0
279       nfl=0                
280       call intout
281 c      call intcartderiv
282 c      call checkintcartgrad
283       call zerograd
284       aincr=1.0D-5
285       write(iout,*) 'Calling CHECK_ECARTINT.'
286       nf=0
287       icall=0
288       call geom_to_var(nvar,x)
289       if (.not.split_ene) then
290         call etotal(energia(0))
291         etot=energia(0)
292         call enerprint(energia(0))
293         call flush(iout)
294         write (iout,*) "enter cartgrad"
295         call flush(iout)
296         call cartgrad
297         write (iout,*) "exit cartgrad"
298         call flush(iout)
299         icall =1
300         do i=1,nres
301           write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
302         enddo
303         do j=1,3
304           grad_s(j,0)=gcart(j,0)
305         enddo
306         do i=1,nres
307           do j=1,3
308             grad_s(j,i)=gcart(j,i)
309             grad_s(j+3,i)=gxcart(j,i)
310           enddo
311         enddo
312       else
313 !- split gradient check
314          write (iout,*) "split_ene not supported"
315 c        call zerograd
316 c        call etotal_long(energia(0))
317 c        call enerprint(energia(0))
318 c        call flush(iout)
319 c        write (iout,*) "enter cartgrad"
320 c        call flush(iout)
321 c        call cartgrad
322 c        write (iout,*) "exit cartgrad"
323 c        call flush(iout)
324 c        icall =1
325 c        write (iout,*) "longrange grad"
326 c        do i=1,nres
327 c          write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
328 c     &    (gxcart(j,i),j=1,3)
329 c        enddo
330 c        do j=1,3
331 c          grad_s(j,0)=gcart(j,0)
332 c        enddo
333 c        do i=1,nres
334 c          do j=1,3
335 c            grad_s(j,i)=gcart(j,i)
336 c            grad_s(j+3,i)=gxcart(j,i)
337 c          enddo
338 c        enddo
339 c        call zerograd
340 c        call etotal_short(energia(0))
341 c        call enerprint(energia(0))
342 c        call flush(iout)
343 c        write (iout,*) "enter cartgrad"
344 c        call flush(iout)
345 c        call cartgrad
346 c        write (iout,*) "exit cartgrad"
347 c        call flush(iout)
348 c        icall =1
349 c        write (iout,*) "shortrange grad"
350 c        do i=1,nres
351 c          write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
352 c     &    (gxcart(j,i),j=1,3)
353 c        enddo
354 c        do j=1,3
355 c          grad_s1(j,0)=gcart(j,0)
356 c        enddo
357 c        do i=1,nres
358 c          do j=1,3
359 c            grad_s1(j,i)=gcart(j,i)
360 c            grad_s1(j+3,i)=gxcart(j,i)
361 c          enddo
362 c        enddo
363       endif
364       write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
365       do i=0,nres
366         do j=1,3
367           xx(j)=c(j,i+nres)
368           ddc(j)=dc(j,i) 
369           ddx(j)=dc(j,i+nres)
370           do k=1,3
371             dcnorm_safe(k)=dc_norm(k,i)
372             dxnorm_safe(k)=dc_norm(k,i+nres)
373           enddo
374         enddo
375         do j=1,3
376           dc(j,i)=ddc(j)+aincr
377           call chainbuild_cart
378 #ifdef MPI
379 c Broadcast the order to compute internal coordinates to the slaves.
380 c          if (nfgtasks.gt.1)
381 c     &      call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
382 #endif
383 c          call int_from_cart1(.false.)
384           if (.not.split_ene) then
385             call etotal(energia1(0))
386             etot1=energia1(0)
387           else
388 !- split gradient
389 c            call etotal_long(energia1(0))
390 c            etot11=energia1(0)
391 c            call etotal_short(energia1(0))
392 c            etot12=energia1(0)
393 c            write (iout,*) "etot11",etot11," etot12",etot12
394           endif
395 !- end split gradient
396 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
397           dc(j,i)=ddc(j)-aincr
398           call chainbuild_cart
399 c          call int_from_cart1(.false.)
400           if (.not.split_ene) then
401             call etotal(energia1(0))
402             etot2=energia1(0)
403             ggg(j)=(etot1-etot2)/(2*aincr)
404           else
405 !- split gradient
406 c            call etotal_long(energia1(0))
407 c            etot21=energia1(0)
408 c           ggg(j)=(etot11-etot21)/(2*aincr)
409 c            call etotal_short(energia1(0))
410 c            etot22=energia1(0)
411 c           ggg1(j)=(etot12-etot22)/(2*aincr)
412 !- end split gradient
413 c            write (iout,*) "etot21",etot21," etot22",etot22
414           endif
415 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
416           dc(j,i)=ddc(j)
417           call chainbuild_cart
418         enddo
419         do j=1,3
420           dc(j,i+nres)=ddx(j)+aincr
421           call chainbuild_cart
422 c          write (iout,*) "i",i," j",j," dxnorm+ and dxnorm"
423 c          write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
424 c          write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
425 c          write (iout,*) "dxnormnorm",dsqrt(
426 c     &  dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
427 c          write (iout,*) "dxnormnormsafe",dsqrt(
428 c     &      dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
429 c          write (iout,*)
430           if (.not.split_ene) then
431             call etotal(energia1(0))
432             etot1=energia1(0)
433           else
434 !- split gradient
435 c            call etotal_long(energia1(0))
436 c            etot11=energia1(0)
437 c            call etotal_short(energia1(0))
438 c            etot12=energia1(0)
439           endif
440 !- end split gradient
441 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
442           dc(j,i+nres)=ddx(j)-aincr
443           call chainbuild_cart
444 c          write (iout,*) "i",i," j",j," dxnorm- and dxnorm"
445 c          write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
446 c          write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
447 c          write (iout,*) 
448 c          write (iout,*) "dxnormnorm",dsqrt(
449 c     &  dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
450 c          write (iout,*) "dxnormnormsafe",dsqrt(
451 c     &      dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
452           if (.not.split_ene) then
453             call etotal(energia1(0))
454             etot2=energia1(0)
455             ggg(j+3)=(etot1-etot2)/(2*aincr)
456           else
457 !- split gradient
458 c            call etotal_long(energia1(0))
459 c            etot21=energia1(0)
460 c           ggg(j+3)=(etot11-etot21)/(2*aincr)
461 c           call etotal_short(energia1(0))
462 c            etot22=energia1(0)
463 c           ggg1(j+3)=(etot12-etot22)/(2*aincr)
464 !- end split gradient
465           endif
466 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
467           dc(j,i+nres)=ddx(j)
468           call chainbuild_cart
469         enddo
470         write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
471      &   i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
472         if (split_ene) then
473           write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
474      &   i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
475      &   k=1,6)
476          write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
477      &   i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
478      &   ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
479         endif
480       enddo
481       return
482       end
483 c-------------------------------------------------------------------------
484       subroutine int_from_cart1(lprn)
485       implicit real*8 (a-h,o-z)
486       include 'DIMENSIONS'
487 #ifdef MPI
488       include 'mpif.h'
489       integer ierror
490 #endif
491       include 'COMMON.IOUNITS'
492       include 'COMMON.VAR'
493       include 'COMMON.CHAIN'
494       include 'COMMON.GEO'
495       include 'COMMON.INTERACT'
496       include 'COMMON.LOCAL'
497       include 'COMMON.NAMES'
498       include 'COMMON.SETUP'
499       include 'COMMON.TIME1'
500       logical lprn 
501       if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
502 #ifdef TIMING
503       time01=MPI_Wtime()
504 #endif
505 #if defined(PARINT) && defined(MPI)
506       do i=iint_start,iint_end
507 #else
508       do i=2,nres
509 #endif
510         dnorm1=dist(i-1,i)
511         dnorm2=dist(i,i+1) 
512         do j=1,3
513           c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1
514      &     +(c(j,i+1)-c(j,i))/dnorm2)
515         enddo
516         be=0.0D0
517         if (i.gt.2) then
518         if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
519         if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
520          tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
521         endif
522         if (itype(i-1).ne.10) then
523          tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
524          omicron(1,i)=alpha(i-2,i-1,i-1+nres)
525          omicron(2,i)=alpha(i-1+nres,i-1,i)
526         endif
527         if (itype(i).ne.10) then
528          tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
529         endif
530         endif
531         omeg(i)=beta(nres+i,i,maxres2,i+1)
532         alph(i)=alpha(nres+i,i,maxres2)
533         theta(i+1)=alpha(i-1,i,i+1)
534         vbld(i)=dist(i-1,i)
535         vbld_inv(i)=1.0d0/vbld(i)
536         vbld(nres+i)=dist(nres+i,i)
537         if (itype(i).ne.10) then
538           vbld_inv(nres+i)=1.0d0/vbld(nres+i)
539         else
540           vbld_inv(nres+i)=0.0d0
541         endif
542       enddo   
543 #if defined(PARINT) && defined(MPI)
544        if (nfgtasks1.gt.1) then
545 cd       write(iout,*) "iint_start",iint_start," iint_count",
546 cd     &   (iint_count(i),i=0,nfgtasks-1)," iint_displ",
547 cd     &   (iint_displ(i),i=0,nfgtasks-1)
548 cd       write (iout,*) "Gather vbld backbone"
549 cd       call flush(iout)
550        time00=MPI_Wtime()
551        call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
552      &   MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
553      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
554 cd       write (iout,*) "Gather vbld_inv"
555 cd       call flush(iout)
556        call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
557      &   MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
558      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
559 cd       write (iout,*) "Gather vbld side chain"
560 cd       call flush(iout)
561        call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
562      &  MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
563      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
564 cd       write (iout,*) "Gather vbld_inv side chain"
565 cd       call flush(iout)
566        call MPI_Allgatherv(vbld_inv(iint_start+nres),
567      &   iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
568      &   iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
569 cd       write (iout,*) "Gather theta"
570 cd       call flush(iout)
571        call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
572      &   MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
573      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
574 cd       write (iout,*) "Gather phi"
575 cd       call flush(iout)
576        call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
577      &   MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
578      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
579 #ifdef CRYST_SC
580 cd       write (iout,*) "Gather alph"
581 cd       call flush(iout)
582        call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
583      &   MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
584      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
585 cd       write (iout,*) "Gather omeg"
586 cd       call flush(iout)
587        call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
588      &   MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
589      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
590 #endif
591        time_gather=time_gather+MPI_Wtime()-time00
592       endif
593 #endif
594       do i=1,nres-1
595         do j=1,3
596           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
597         enddo
598       enddo
599       do i=2,nres-1
600         do j=1,3
601           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
602         enddo
603       enddo
604       if (lprn) then
605       do i=2,nres
606        write (iout,1212) restyp(itype(i)),i,vbld(i),
607      &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
608      &rad2deg*alph(i),rad2deg*omeg(i)
609       enddo
610       endif
611  1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
612 #ifdef TIMING
613       time_intfcart=time_intfcart+MPI_Wtime()-time01
614 #endif
615       return
616       end
617 c----------------------------------------------------------------------------
618       subroutine check_eint
619 C Check the gradient of energy in internal coordinates.
620       implicit real*8 (a-h,o-z)
621       include 'DIMENSIONS'
622       include 'COMMON.CHAIN'
623       include 'COMMON.DERIV'
624       include 'COMMON.IOUNITS'
625       include 'COMMON.VAR'
626       include 'COMMON.GEO'
627       common /srutu/ icall
628       dimension x(maxvar),gana(maxvar),gg(maxvar)
629       integer uiparm(1)
630       double precision urparm(1)
631       double precision energia(0:n_ene),energia1(0:n_ene),
632      &  energia2(0:n_ene)
633       character*6 key
634       external fdum
635       call zerograd
636       aincr=1.0D-7
637       print '(a)','Calling CHECK_INT.'
638       nf=0
639       nfl=0
640       icg=1
641       call geom_to_var(nvar,x)
642       call var_to_geom(nvar,x)
643       call chainbuild
644       icall=1
645       print *,'ICG=',ICG
646       call etotal(energia(0))
647       etot = energia(0)
648       call enerprint(energia(0))
649       print *,'ICG=',ICG
650 #ifdef MPL
651       if (MyID.ne.BossID) then
652         call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
653         nf=x(nvar+1)
654         nfl=x(nvar+2)
655         icg=x(nvar+3)
656       endif
657 #endif
658       nf=1
659       nfl=3
660 cd    write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
661       call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
662 cd    write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
663       icall=1
664       do i=1,nvar
665         xi=x(i)
666         x(i)=xi-0.5D0*aincr
667         call var_to_geom(nvar,x)
668         call chainbuild
669         call etotal(energia1(0))
670         etot1=energia1(0)
671         x(i)=xi+0.5D0*aincr
672         call var_to_geom(nvar,x)
673         call chainbuild
674         call etotal(energia2(0))
675         etot2=energia2(0)
676         gg(i)=(etot2-etot1)/aincr
677         write (iout,*) i,etot1,etot2
678         x(i)=xi
679       enddo
680       write (iout,'(/2a)')' Variable        Numerical       Analytical',
681      &    '     RelDiff*100% '
682       do i=1,nvar
683         if (i.le.nphi) then
684           ii=i
685           key = ' phi'
686         else if (i.le.nphi+ntheta) then
687           ii=i-nphi
688           key=' theta'
689         else if (i.le.nphi+ntheta+nside) then
690            ii=i-(nphi+ntheta)
691            key=' alpha'
692         else 
693            ii=i-(nphi+ntheta+nside)
694            key=' omega'
695         endif
696         write (iout,'(i3,a,i3,3(1pd16.6))') 
697      & i,key,ii,gg(i),gana(i),
698      & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)
699       enddo
700       return
701       end