wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / source / unres / src_MD-M / 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 c      r_cut=2.0d0
276 c      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-4
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         call zerograd
315         call etotal_long(energia(0))
316         call enerprint(energia(0))
317         call flush(iout)
318         write (iout,*) "enter cartgrad"
319         call flush(iout)
320         call cartgrad
321         write (iout,*) "exit cartgrad"
322         call flush(iout)
323         icall =1
324         write (iout,*) "longrange grad"
325         do i=1,nres
326           write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
327      &    (gxcart(j,i),j=1,3)
328         enddo
329         do j=1,3
330           grad_s(j,0)=gcart(j,0)
331         enddo
332         do i=1,nres
333           do j=1,3
334             grad_s(j,i)=gcart(j,i)
335             grad_s(j+3,i)=gxcart(j,i)
336           enddo
337         enddo
338         call zerograd
339         call etotal_short(energia(0))
340         call enerprint(energia(0))
341         call flush(iout)
342         write (iout,*) "enter cartgrad"
343         call flush(iout)
344         call cartgrad
345         write (iout,*) "exit cartgrad"
346         call flush(iout)
347         icall =1
348         write (iout,*) "shortrange grad"
349         do i=1,nres
350           write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
351      &    (gxcart(j,i),j=1,3)
352         enddo
353         do j=1,3
354           grad_s1(j,0)=gcart(j,0)
355         enddo
356         do i=1,nres
357           do j=1,3
358             grad_s1(j,i)=gcart(j,i)
359             grad_s1(j+3,i)=gxcart(j,i)
360           enddo
361         enddo
362       endif
363       write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
364       do i=0,nres
365         do j=1,3
366           xx(j)=c(j,i+nres)
367           ddc(j)=dc(j,i) 
368           ddx(j)=dc(j,i+nres)
369           do k=1,3
370             dcnorm_safe(k)=dc_norm(k,i)
371             dxnorm_safe(k)=dc_norm(k,i+nres)
372           enddo
373         enddo
374         do j=1,3
375           dc(j,i)=ddc(j)+aincr
376           call chainbuild_cart
377 #ifdef MPI
378 c Broadcast the order to compute internal coordinates to the slaves.
379 c          if (nfgtasks.gt.1)
380 c     &      call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
381 #endif
382 c          call int_from_cart1(.false.)
383           if (.not.split_ene) then
384             call etotal(energia1(0))
385             etot1=energia1(0)
386           else
387 !- split gradient
388             call etotal_long(energia1(0))
389             etot11=energia1(0)
390             call etotal_short(energia1(0))
391             etot12=energia1(0)
392 c            write (iout,*) "etot11",etot11," etot12",etot12
393           endif
394 !- end split gradient
395 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
396           dc(j,i)=ddc(j)-aincr
397           call chainbuild_cart
398 C          print *,c(j,i)
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             call etotal_long(energia1(0))
407             etot21=energia1(0)
408             ggg(j)=(etot11-etot21)/(2*aincr)
409             call etotal_short(energia1(0))
410             etot22=energia1(0)
411             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             call etotal_long(energia1(0))
436             etot11=energia1(0)
437             call etotal_short(energia1(0))
438             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             call etotal_long(energia1(0))
459             etot21=energia1(0)
460             ggg(j+3)=(etot11-etot21)/(2*aincr)
461             call etotal_short(energia1(0))
462             etot22=energia1(0)
463             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