ctest the first static disulfide test ene/min cart
[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       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-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          call int_from_cart1(.false.)
399           if (.not.split_ene) then
400             call etotal(energia1(0))
401             etot2=energia1(0)
402             ggg(j)=(etot1-etot2)/(2*aincr)
403           else
404 !- split gradient
405             call etotal_long(energia1(0))
406             etot21=energia1(0)
407             ggg(j)=(etot11-etot21)/(2*aincr)
408             call etotal_short(energia1(0))
409             etot22=energia1(0)
410             ggg1(j)=(etot12-etot22)/(2*aincr)
411 !- end split gradient
412 c            write (iout,*) "etot21",etot21," etot22",etot22
413           endif
414 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
415           dc(j,i)=ddc(j)
416           call chainbuild_cart
417         enddo
418         do j=1,3
419           dc(j,i+nres)=ddx(j)+aincr
420           call chainbuild_cart
421 c          write (iout,*) "i",i," j",j," dxnorm+ and dxnorm"
422 c          write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
423 c          write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
424 c          write (iout,*) "dxnormnorm",dsqrt(
425 c     &  dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
426 c          write (iout,*) "dxnormnormsafe",dsqrt(
427 c     &      dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
428 c          write (iout,*)
429           if (.not.split_ene) then
430             call etotal(energia1(0))
431             etot1=energia1(0)
432           else
433 !- split gradient
434             call etotal_long(energia1(0))
435             etot11=energia1(0)
436             call etotal_short(energia1(0))
437             etot12=energia1(0)
438           endif
439 !- end split gradient
440 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
441           dc(j,i+nres)=ddx(j)-aincr
442           call chainbuild_cart
443 c          write (iout,*) "i",i," j",j," dxnorm- and dxnorm"
444 c          write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
445 c          write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
446 c          write (iout,*) 
447 c          write (iout,*) "dxnormnorm",dsqrt(
448 c     &  dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
449 c          write (iout,*) "dxnormnormsafe",dsqrt(
450 c     &      dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
451           if (.not.split_ene) then
452             call etotal(energia1(0))
453             etot2=energia1(0)
454             ggg(j+3)=(etot1-etot2)/(2*aincr)
455           else
456 !- split gradient
457             call etotal_long(energia1(0))
458             etot21=energia1(0)
459             ggg(j+3)=(etot11-etot21)/(2*aincr)
460             call etotal_short(energia1(0))
461             etot22=energia1(0)
462             ggg1(j+3)=(etot12-etot22)/(2*aincr)
463 !- end split gradient
464           endif
465 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
466           dc(j,i+nres)=ddx(j)
467           call chainbuild_cart
468         enddo
469         write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
470      &   i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
471         if (split_ene) then
472           write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
473      &   i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
474      &   k=1,6)
475          write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
476      &   i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
477      &   ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
478         endif
479       enddo
480       return
481       end
482 c-------------------------------------------------------------------------
483       subroutine int_from_cart1(lprn)
484       implicit real*8 (a-h,o-z)
485       include 'DIMENSIONS'
486 #ifdef MPI
487       include 'mpif.h'
488       integer ierror
489 #endif
490       include 'COMMON.IOUNITS'
491       include 'COMMON.VAR'
492       include 'COMMON.CHAIN'
493       include 'COMMON.GEO'
494       include 'COMMON.INTERACT'
495       include 'COMMON.LOCAL'
496       include 'COMMON.NAMES'
497       include 'COMMON.SETUP'
498       include 'COMMON.TIME1'
499       logical lprn 
500       if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
501 #ifdef TIMING
502       time01=MPI_Wtime()
503 #endif
504 #if defined(PARINT) && defined(MPI)
505       do i=iint_start,iint_end
506 #else
507       do i=2,nres
508 #endif
509         dnorm1=dist(i-1,i)
510         dnorm2=dist(i,i+1) 
511         do j=1,3
512           c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1
513      &     +(c(j,i+1)-c(j,i))/dnorm2)
514         enddo
515         be=0.0D0
516         if (i.gt.2) then
517         if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
518         if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
519          tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
520         endif
521         if (itype(i-1).ne.10) then
522          tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
523          omicron(1,i)=alpha(i-2,i-1,i-1+nres)
524          omicron(2,i)=alpha(i-1+nres,i-1,i)
525         endif
526         if (itype(i).ne.10) then
527          tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
528         endif
529         endif
530         omeg(i)=beta(nres+i,i,maxres2,i+1)
531         alph(i)=alpha(nres+i,i,maxres2)
532         theta(i+1)=alpha(i-1,i,i+1)
533         vbld(i)=dist(i-1,i)
534         vbld_inv(i)=1.0d0/vbld(i)
535         vbld(nres+i)=dist(nres+i,i)
536         if (itype(i).ne.10) then
537           vbld_inv(nres+i)=1.0d0/vbld(nres+i)
538         else
539           vbld_inv(nres+i)=0.0d0
540         endif
541       enddo   
542 #if defined(PARINT) && defined(MPI)
543        if (nfgtasks1.gt.1) then
544 cd       write(iout,*) "iint_start",iint_start," iint_count",
545 cd     &   (iint_count(i),i=0,nfgtasks-1)," iint_displ",
546 cd     &   (iint_displ(i),i=0,nfgtasks-1)
547 cd       write (iout,*) "Gather vbld backbone"
548 cd       call flush(iout)
549        time00=MPI_Wtime()
550        call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
551      &   MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
552      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
553 cd       write (iout,*) "Gather vbld_inv"
554 cd       call flush(iout)
555        call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
556      &   MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
557      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
558 cd       write (iout,*) "Gather vbld side chain"
559 cd       call flush(iout)
560        call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
561      &  MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
562      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
563 cd       write (iout,*) "Gather vbld_inv side chain"
564 cd       call flush(iout)
565        call MPI_Allgatherv(vbld_inv(iint_start+nres),
566      &   iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
567      &   iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
568 cd       write (iout,*) "Gather theta"
569 cd       call flush(iout)
570        call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
571      &   MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
572      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
573 cd       write (iout,*) "Gather phi"
574 cd       call flush(iout)
575        call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
576      &   MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
577      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
578 #ifdef CRYST_SC
579 cd       write (iout,*) "Gather alph"
580 cd       call flush(iout)
581        call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
582      &   MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
583      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
584 cd       write (iout,*) "Gather omeg"
585 cd       call flush(iout)
586        call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
587      &   MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
588      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
589 #endif
590        time_gather=time_gather+MPI_Wtime()-time00
591       endif
592 #endif
593       do i=1,nres-1
594         do j=1,3
595           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
596         enddo
597       enddo
598       do i=2,nres-1
599         do j=1,3
600           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
601         enddo
602       enddo
603       if (lprn) then
604       do i=2,nres
605        write (iout,1212) restyp(itype(i)),i,vbld(i),
606      &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
607      &rad2deg*alph(i),rad2deg*omeg(i)
608       enddo
609       endif
610  1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
611 #ifdef TIMING
612       time_intfcart=time_intfcart+MPI_Wtime()-time01
613 #endif
614       return
615       end
616 c----------------------------------------------------------------------------
617       subroutine check_eint
618 C Check the gradient of energy in internal coordinates.
619       implicit real*8 (a-h,o-z)
620       include 'DIMENSIONS'
621       include 'COMMON.CHAIN'
622       include 'COMMON.DERIV'
623       include 'COMMON.IOUNITS'
624       include 'COMMON.VAR'
625       include 'COMMON.GEO'
626       common /srutu/ icall
627       dimension x(maxvar),gana(maxvar),gg(maxvar)
628       integer uiparm(1)
629       double precision urparm(1)
630       double precision energia(0:n_ene),energia1(0:n_ene),
631      &  energia2(0:n_ene)
632       character*6 key
633       external fdum
634       call zerograd
635       aincr=1.0D-7
636       print '(a)','Calling CHECK_INT.'
637       nf=0
638       nfl=0
639       icg=1
640       call geom_to_var(nvar,x)
641       call var_to_geom(nvar,x)
642       call chainbuild
643       icall=1
644       print *,'ICG=',ICG
645       call etotal(energia(0))
646       etot = energia(0)
647       call enerprint(energia(0))
648       print *,'ICG=',ICG
649 #ifdef MPL
650       if (MyID.ne.BossID) then
651         call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
652         nf=x(nvar+1)
653         nfl=x(nvar+2)
654         icg=x(nvar+3)
655       endif
656 #endif
657       nf=1
658       nfl=3
659 cd    write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
660       call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
661 cd    write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
662       icall=1
663       do i=1,nvar
664         xi=x(i)
665         x(i)=xi-0.5D0*aincr
666         call var_to_geom(nvar,x)
667         call chainbuild
668         call etotal(energia1(0))
669         etot1=energia1(0)
670         x(i)=xi+0.5D0*aincr
671         call var_to_geom(nvar,x)
672         call chainbuild
673         call etotal(energia2(0))
674         etot2=energia2(0)
675         gg(i)=(etot2-etot1)/aincr
676         write (iout,*) i,etot1,etot2
677         x(i)=xi
678       enddo
679       write (iout,'(/2a)')' Variable        Numerical       Analytical',
680      &    '     RelDiff*100% '
681       do i=1,nvar
682         if (i.le.nphi) then
683           ii=i
684           key = ' phi'
685         else if (i.le.nphi+ntheta) then
686           ii=i-nphi
687           key=' theta'
688         else if (i.le.nphi+ntheta+nside) then
689            ii=i-(nphi+ntheta)
690            key=' alpha'
691         else 
692            ii=i-(nphi+ntheta+nside)
693            key=' omega'
694         endif
695         write (iout,'(i3,a,i3,3(1pd16.6))') 
696      & i,key,ii,gg(i),gana(i),
697      & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)
698       enddo
699       return
700       end