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