Merge branch 'bartek' of mmka.chem.univ.gda.pl:unres into bartek
[unres.git] / source / unres / src_MD / intcartderiv.F
1       subroutine intcartderiv
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4 #ifdef MPI
5       include 'mpif.h'
6 #endif
7       include 'COMMON.SETUP'
8       include 'COMMON.CHAIN' 
9       include 'COMMON.VAR'
10       include 'COMMON.GEO'
11       include 'COMMON.INTERACT'
12       include 'COMMON.DERIV'
13       include 'COMMON.IOUNITS'
14       include 'COMMON.LOCAL'
15       include 'COMMON.SCCOR'
16       double precision dcostheta(3,2,maxres),
17      & dcosphi(3,3,maxres),dsinphi(3,3,maxres),
18      & dcosalpha(3,3,maxres),dcosomega(3,3,maxres),
19      & dsinomega(3,3,maxres),vo1(3),vo2(3),vo3(3),
20      & dummy(3),vp1(3),vp2(3),vp3(3),vpp1(3),n(3)
21        
22 #if defined(MPI) && defined(PARINTDER)
23       if (nfgtasks.gt.1 .and. me.eq.king) 
24      &  call MPI_Bcast(8,1,MPI_INTEGER,king,FG_COMM,IERROR)
25 #endif
26       pi4 = 0.5d0*pipol
27       pi34 = 3*pi4
28       
29 c      write (iout,*) "iphi1_start",iphi1_start," iphi1_end",iphi1_end      
30 c Derivatives of theta's
31 #if defined(MPI) && defined(PARINTDER)
32 c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
33       do i=max0(ithet_start-1,3),ithet_end
34 #else
35       do i=3,nres
36 #endif
37         cost=dcos(theta(i))
38         sint=sqrt(1-cost*cost)
39         do j=1,3
40           dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/
41      &    vbld(i-1)
42           dtheta(j,1,i)=-1/sint*dcostheta(j,1,i)     
43           dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/
44      &    vbld(i)
45           dtheta(j,2,i)=-1/sint*dcostheta(j,2,i)     
46         enddo
47       enddo
48
49 #if defined(MPI) && defined(PARINTDER)
50 c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
51       do i=max0(ithet_start-1,3),ithet_end
52 #else
53       do i=3,nres
54 #endif
55       if ((itype(i-1).ne.10).and.(itype(i-1).ne.21)) then
56         cost1=dcos(omicron(1,i))
57         sint1=sqrt(1-cost1*cost1)
58         cost2=dcos(omicron(2,i))
59         sint2=sqrt(1-cost2*cost2)
60         do j=1,3
61 CC Calculate derivative over first omicron (Cai-2,Cai-1,SCi-1) 
62           dcosomicron(j,1,1,i)=-(dc_norm(j,i-1+nres)+
63      &    cost1*dc_norm(j,i-2))/
64      &    vbld(i-1)
65           domicron(j,1,1,i)=-1/sint1*dcosomicron(j,1,1,i)     
66           dcosomicron(j,1,2,i)=-(dc_norm(j,i-2)
67      &    +cost1*(dc_norm(j,i-1+nres)))/
68      &    vbld(i-1+nres)
69           domicron(j,1,2,i)=-1/sint1*dcosomicron(j,1,2,i)  
70 CC Calculate derivative over second omicron Sci-1,Cai-1 Cai
71 CC Looks messy but better than if in loop
72           dcosomicron(j,2,1,i)=-(-dc_norm(j,i-1+nres)
73      &    +cost2*dc_norm(j,i-1))/
74      &    vbld(i)
75           domicron(j,2,1,i)=-1/sint2*dcosomicron(j,2,1,i)
76           dcosomicron(j,2,2,i)=-(dc_norm(j,i-1)
77      &     +cost2*(-dc_norm(j,i-1+nres)))/
78      &    vbld(i-1+nres)
79 c          write(iout,*) "vbld", i,itype(i),vbld(i-1+nres)
80           domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i)   
81         enddo
82        endif
83       enddo
84
85
86       
87 c Derivatives of phi:
88 c If phi is 0 or 180 degrees, then the formulas 
89 c have to be derived by power series expansion of the
90 c conventional formulas around 0 and 180.
91 #ifdef PARINTDER
92       do i=iphi1_start,iphi1_end
93 #else
94       do i=4,nres      
95 #endif
96 c the conventional case
97         sint=dsin(theta(i))
98         sint1=dsin(theta(i-1))
99         sing=dsin(phi(i))
100         cost=dcos(theta(i))
101         cost1=dcos(theta(i-1))
102         cosg=dcos(phi(i))
103         scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1))
104         fac0=1.0d0/(sint1*sint)
105         fac1=cost*fac0
106         fac2=cost1*fac0
107         fac3=cosg*cost1/(sint1*sint1)
108         fac4=cosg*cost/(sint*sint)
109 c    Obtaining the gamma derivatives from sine derivative                                
110        if (phi(i).gt.-pi4.and.phi(i).le.pi4.or.
111      &     phi(i).gt.pi34.and.phi(i).le.pi.or.
112      &     phi(i).gt.-pi.and.phi(i).le.-pi34) then
113          call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
114          call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2)
115          call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3) 
116          do j=1,3
117             ctgt=cost/sint
118             ctgt1=cost1/sint1
119             cosg_inv=1.0d0/cosg
120             dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
121      &        -(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
122             dphi(j,1,i)=cosg_inv*dsinphi(j,1,i)
123             dsinphi(j,2,i)=
124      &        -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*dtheta(j,1,i))
125      &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
126             dphi(j,2,i)=cosg_inv*dsinphi(j,2,i)
127 c Bug fixed 3/24/05 (AL)
128             dsinphi(j,3,i)=-sing*ctgt*dtheta(j,2,i)
129      &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
130 c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
131             dphi(j,3,i)=cosg_inv*dsinphi(j,3,i)
132          enddo                                              
133 c   Obtaining the gamma derivatives from cosine derivative
134         else
135            do j=1,3
136            dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
137      &     dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp*
138      &     dc_norm(j,i-3))/vbld(i-2)
139            dphi(j,1,i)=-1/sing*dcosphi(j,1,i)       
140            dcosphi(j,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
141      &     dcostheta(j,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
142      &     dcostheta(j,1,i)
143            dphi(j,2,i)=-1/sing*dcosphi(j,2,i)      
144            dcosphi(j,3,i)=fac2*dcostheta(j,2,i)+fac4*
145      &     dcostheta(j,2,i)-fac0*(dc_norm(j,i-3)-scalp*
146      &     dc_norm(j,i-1))/vbld(i)
147            dphi(j,3,i)=-1/sing*dcosphi(j,3,i)       
148          enddo
149         endif                                                                                            
150       enddo
151
152 Calculate derivative of Tauangle
153 #ifdef PARINTDER
154       do i=itau_start,itau_end
155 #else
156       do i=3,nres
157 #endif
158        if ((itype(i-2).eq.21).or.(itype(i-2).eq.10)) cycle
159 cc dtauangle(j,intertyp,dervityp,residue number)
160 cc INTERTYP=1 SC...Ca...Ca..Ca
161 c the conventional case
162         sint=dsin(theta(i))
163         sint1=dsin(omicron(2,i-1))
164         sing=dsin(tauangle(1,i))
165         cost=dcos(theta(i))
166         cost1=dcos(omicron(2,i-1))
167         cosg=dcos(tauangle(1,i))
168         do j=1,3
169         dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
170 cc       write(iout,*) dc_norm2(j,i-2+nres),"dcnorm"
171         enddo
172         scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1))
173         fac0=1.0d0/(sint1*sint)
174         fac1=cost*fac0
175         fac2=cost1*fac0
176         fac3=cosg*cost1/(sint1*sint1)
177         fac4=cosg*cost/(sint*sint)
178 cc         write(iout,*) "faki",fac0,fac1,fac2,fac3,fac4
179 c    Obtaining the gamma derivatives from sine derivative                                
180        if (tauangle(1,i).gt.-pi4.and.tauangle(1,i).le.pi4.or.
181      &     tauangle(1,i).gt.pi34.and.tauangle(1,i).le.pi.or.
182      &     tauangle(1,i).gt.-pi.and.tauangle(1,i).le.-pi34) then
183          call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
184          call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1),vp2)
185          call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
186         do j=1,3
187             ctgt=cost/sint
188             ctgt1=cost1/sint1
189             cosg_inv=1.0d0/cosg
190             dsintau(j,1,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
191      &-(fac0*vp1(j)+sing*(dc_norm2(j,i-2+nres)))
192      & *vbld_inv(i-2+nres)
193             dtauangle(j,1,1,i)=cosg_inv*dsintau(j,1,1,i)
194             dsintau(j,1,2,i)=
195      &        -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*dtheta(j,1,i))
196      &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
197 c            write(iout,*) "dsintau", dsintau(j,1,2,i)
198             dtauangle(j,1,2,i)=cosg_inv*dsintau(j,1,2,i)
199 c Bug fixed 3/24/05 (AL)
200             dsintau(j,1,3,i)=-sing*ctgt*dtheta(j,2,i)
201      &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
202 c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
203             dtauangle(j,1,3,i)=cosg_inv*dsintau(j,1,3,i)
204          enddo                                          
205 c   Obtaining the gamma derivatives from cosine derivative
206         else
207            do j=1,3
208            dcostau(j,1,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
209      &     dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1)-scalp*
210      &     (dc_norm2(j,i-2+nres)))/vbld(i-2+nres)
211            dtauangle(j,1,1,i)=-1/sing*dcostau(j,1,1,i)
212            dcostau(j,1,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
213      &     dcostheta(j,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
214      &     dcostheta(j,1,i)
215            dtauangle(j,1,2,i)=-1/sing*dcostau(j,1,2,i)
216            dcostau(j,1,3,i)=fac2*dcostheta(j,2,i)+fac4*
217      &     dcostheta(j,2,i)-fac0*(-dc_norm(j,i-2+nres)-scalp*
218      &     dc_norm(j,i-1))/vbld(i)
219            dtauangle(j,1,3,i)=-1/sing*dcostau(j,1,3,i)
220 c         write (iout,*) "else",i
221          enddo
222         endif
223 c        do k=1,3                 
224 c        write(iout,*) "tu",i,k,(dtauangle(j,1,k,i),j=1,3)        
225 c        enddo                
226       enddo
227 CC Second case Ca...Ca...Ca...SC
228 #ifdef PARINTDER
229       do i=itau_start,itau_end
230 #else
231       do i=4,nres
232 #endif
233        if ((itype(i-1).eq.21).or.(itype(i-1).eq.10)) cycle
234 c the conventional case
235         sint=dsin(omicron(1,i))
236         sint1=dsin(theta(i-1))
237         sing=dsin(tauangle(2,i))
238         cost=dcos(omicron(1,i))
239         cost1=dcos(theta(i-1))
240         cosg=dcos(tauangle(2,i))
241 c        do j=1,3
242 c        dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
243 c        enddo
244         scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1+nres))
245         fac0=1.0d0/(sint1*sint)
246         fac1=cost*fac0
247         fac2=cost1*fac0
248         fac3=cosg*cost1/(sint1*sint1)
249         fac4=cosg*cost/(sint*sint)
250 c    Obtaining the gamma derivatives from sine derivative                                
251        if (tauangle(2,i).gt.-pi4.and.tauangle(2,i).le.pi4.or.
252      &     tauangle(2,i).gt.pi34.and.tauangle(2,i).le.pi.or.
253      &     tauangle(2,i).gt.-pi.and.tauangle(2,i).le.-pi34) then
254          call vecpr(dc_norm2(1,i-1+nres),dc_norm(1,i-2),vp1)
255          call vecpr(dc_norm(1,i-3),dc_norm(1,i-1+nres),vp2)
256          call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
257         do j=1,3
258             ctgt=cost/sint
259             ctgt1=cost1/sint1
260             cosg_inv=1.0d0/cosg
261             dsintau(j,2,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
262      &        +(fac0*vp1(j)-sing*dc_norm(j,i-3))*vbld_inv(i-2)
263 c       write(iout,*) i,j,dsintau(j,2,1,i),sing*ctgt1*dtheta(j,1,i-1),
264 c     &fac0*vp1(j),sing*dc_norm(j,i-3),vbld_inv(i-2),"dsintau(2,1)"
265             dtauangle(j,2,1,i)=cosg_inv*dsintau(j,2,1,i)
266             dsintau(j,2,2,i)=
267      &        -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*domicron(j,1,1,i))
268      &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
269 c            write(iout,*) "sprawdzenie",i,j,sing*ctgt1*dtheta(j,2,i-1),
270 c     & sing*ctgt*domicron(j,1,2,i),
271 c     & (fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
272             dtauangle(j,2,2,i)=cosg_inv*dsintau(j,2,2,i)
273 c Bug fixed 3/24/05 (AL)
274             dsintau(j,2,3,i)=-sing*ctgt*domicron(j,1,2,i)
275      &       +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))*vbld_inv(i-1+nres)
276 c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
277             dtauangle(j,2,3,i)=cosg_inv*dsintau(j,2,3,i)
278          enddo                                          
279 c   Obtaining the gamma derivatives from cosine derivative
280         else
281            do j=1,3
282            dcostau(j,2,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
283      &     dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
284      &     dc_norm(j,i-3))/vbld(i-2)
285            dtauangle(j,2,1,i)=-1/sing*dcostau(j,2,1,i)
286            dcostau(j,2,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
287      &     dcosomicron(j,1,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
288      &     dcosomicron(j,1,1,i)
289            dtauangle(j,2,2,i)=-1/sing*dcostau(j,2,2,i)
290            dcostau(j,2,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
291      &     dcosomicron(j,1,2,i)-fac0*(dc_norm(j,i-3)-scalp*
292      &     dc_norm(j,i-1+nres))/vbld(i-1+nres)
293            dtauangle(j,2,3,i)=-1/sing*dcostau(j,2,3,i)
294 c        write(iout,*) i,j,"else", dtauangle(j,2,3,i) 
295          enddo
296         endif                                                                                            
297       enddo
298
299
300 CCC third case SC...Ca...Ca...SC
301 #ifdef PARINTDER
302
303       do i=itau_start,itau_end
304 #else
305       do i=3,nres
306 #endif
307 c the conventional case
308       if ((itype(i-1).eq.21).or.(itype(i-1).eq.10).or.
309      &(itype(i-2).eq.21).or.(itype(i-2).eq.10)) cycle
310         sint=dsin(omicron(1,i))
311         sint1=dsin(omicron(2,i-1))
312         sing=dsin(tauangle(3,i))
313         cost=dcos(omicron(1,i))
314         cost1=dcos(omicron(2,i-1))
315         cosg=dcos(tauangle(3,i))
316         do j=1,3
317         dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
318 c        dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
319         enddo
320         scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres))
321         fac0=1.0d0/(sint1*sint)
322         fac1=cost*fac0
323         fac2=cost1*fac0
324         fac3=cosg*cost1/(sint1*sint1)
325         fac4=cosg*cost/(sint*sint)
326 c    Obtaining the gamma derivatives from sine derivative                                
327        if (tauangle(3,i).gt.-pi4.and.tauangle(3,i).le.pi4.or.
328      &     tauangle(3,i).gt.pi34.and.tauangle(3,i).le.pi.or.
329      &     tauangle(3,i).gt.-pi.and.tauangle(3,i).le.-pi34) then
330          call vecpr(dc_norm(1,i-1+nres),dc_norm(1,i-2),vp1)
331          call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres),vp2)
332          call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
333         do j=1,3
334             ctgt=cost/sint
335             ctgt1=cost1/sint1
336             cosg_inv=1.0d0/cosg
337             dsintau(j,3,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
338      &        -(fac0*vp1(j)-sing*dc_norm(j,i-2+nres))
339      &        *vbld_inv(i-2+nres)
340             dtauangle(j,3,1,i)=cosg_inv*dsintau(j,3,1,i)
341             dsintau(j,3,2,i)=
342      &        -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*domicron(j,1,1,i))
343      &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
344             dtauangle(j,3,2,i)=cosg_inv*dsintau(j,3,2,i)
345 c Bug fixed 3/24/05 (AL)
346             dsintau(j,3,3,i)=-sing*ctgt*domicron(j,1,2,i)
347      &        +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))
348      &        *vbld_inv(i-1+nres)
349 c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
350             dtauangle(j,3,3,i)=cosg_inv*dsintau(j,3,3,i)
351          enddo                                          
352 c   Obtaining the gamma derivatives from cosine derivative
353         else
354            do j=1,3
355            dcostau(j,3,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
356      &     dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
357      &     dc_norm2(j,i-2+nres))/vbld(i-2+nres)
358            dtauangle(j,3,1,i)=-1/sing*dcostau(j,3,1,i)
359            dcostau(j,3,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
360      &     dcosomicron(j,1,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
361      &     dcosomicron(j,1,1,i)
362            dtauangle(j,3,2,i)=-1/sing*dcostau(j,3,2,i)
363            dcostau(j,3,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
364      &     dcosomicron(j,1,2,i)-fac0*(dc_norm2(j,i-2+nres)-scalp*
365      &     dc_norm(j,i-1+nres))/vbld(i-1+nres)
366            dtauangle(j,3,3,i)=-1/sing*dcostau(j,3,3,i)
367 c          write(iout,*) "else",i 
368          enddo
369         endif                                                                                            
370       enddo
371 #ifdef CRYST_SC
372 c   Derivatives of side-chain angles alpha and omega
373 #if defined(MPI) && defined(PARINTDER)
374         do i=ibond_start,ibond_end
375 #else
376         do i=2,nres-1           
377 #endif
378           if(itype(i).ne.10) then         
379              fac5=1.0d0/dsqrt(2*(1+dcos(theta(i+1))))
380              fac6=fac5/vbld(i)
381              fac7=fac5*fac5
382              fac8=fac5/vbld(i+1)     
383              fac9=fac5/vbld(i+nres)                  
384              scala1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
385              scala2=scalar(dc_norm(1,i),dc_norm(1,i+nres))
386              cosa=dsqrt(0.5d0/(1.0d0+dcos(theta(i+1))))*(
387      &       scalar(dC_norm(1,i),dC_norm(1,i+nres))
388      &       -scalar(dC_norm(1,i-1),dC_norm(1,i+nres)))
389              sina=sqrt(1-cosa*cosa)
390              sino=dsin(omeg(i))                                                                                              
391              do j=1,3     
392                 dcosalpha(j,1,i)=fac6*(scala1*dc_norm(j,i-1)-
393      &          dc_norm(j,i+nres))-cosa*fac7*dcostheta(j,1,i+1)
394                 dalpha(j,1,i)=-1/sina*dcosalpha(j,1,i)
395                 dcosalpha(j,2,i)=fac8*(dc_norm(j,i+nres)-
396      &          scala2*dc_norm(j,i))-cosa*fac7*dcostheta(j,2,i+1)
397                 dalpha(j,2,i)=-1/sina*dcosalpha(j,2,i)
398                 dcosalpha(j,3,i)=(fac9*(dc_norm(j,i)-
399      &          dc_norm(j,i-1))-(cosa*dc_norm(j,i+nres))/
400      &          vbld(i+nres))
401                 dalpha(j,3,i)=-1/sina*dcosalpha(j,3,i)
402             enddo
403 c obtaining the derivatives of omega from sines     
404             if(omeg(i).gt.-pi4.and.omeg(i).le.pi4.or.
405      &         omeg(i).gt.pi34.and.omeg(i).le.pi.or.
406      &         omeg(i).gt.-pi.and.omeg(i).le.-pi34) then
407                fac15=dcos(theta(i+1))/(dsin(theta(i+1))*
408      &         dsin(theta(i+1)))
409                fac16=dcos(alph(i))/(dsin(alph(i))*dsin(alph(i)))
410                fac17=1.0d0/(dsin(theta(i+1))*dsin(alph(i)))             
411                call vecpr(dc_norm(1,i+nres),dc_norm(1,i),vo1)
412                call vecpr(dc_norm(1,i+nres),dc_norm(1,i-1),vo2)
413                call vecpr(dc_norm(1,i),dc_norm(1,i-1),vo3)
414                coso_inv=1.0d0/dcos(omeg(i))                            
415                do j=1,3
416                  dsinomega(j,1,i)=sino*(fac15*dcostheta(j,1,i+1)
417      &           +fac16*dcosalpha(j,1,i))-fac17/vbld(i)*vo1(j)-(
418      &           sino*dc_norm(j,i-1))/vbld(i)
419                  domega(j,1,i)=coso_inv*dsinomega(j,1,i)
420                  dsinomega(j,2,i)=sino*(fac15*dcostheta(j,2,i+1)
421      &           +fac16*dcosalpha(j,2,i))+fac17/vbld(i+1)*vo2(j)
422      &           -sino*dc_norm(j,i)/vbld(i+1)
423                  domega(j,2,i)=coso_inv*dsinomega(j,2,i)                                                       
424                  dsinomega(j,3,i)=sino*fac16*dcosalpha(j,3,i)-
425      &           fac17/vbld(i+nres)*vo3(j)-sino*dc_norm(j,i+nres)/
426      &           vbld(i+nres)
427                  domega(j,3,i)=coso_inv*dsinomega(j,3,i)
428               enddo                              
429            else
430 c   obtaining the derivatives of omega from cosines
431              fac10=sqrt(0.5d0*(1-dcos(theta(i+1))))
432              fac11=sqrt(0.5d0*(1+dcos(theta(i+1))))
433              fac12=fac10*sina
434              fac13=fac12*fac12
435              fac14=sina*sina
436              do j=1,3                                    
437                 dcosomega(j,1,i)=(-(0.25d0*cosa/fac11*
438      &          dcostheta(j,1,i+1)+fac11*dcosalpha(j,1,i))*fac12+
439      &          (0.25d0/fac10*sina*dcostheta(j,1,i+1)+cosa/sina*
440      &          fac10*dcosalpha(j,1,i))*(scala2-fac11*cosa))/fac13
441                 domega(j,1,i)=-1/sino*dcosomega(j,1,i)
442                 dcosomega(j,2,i)=(((dc_norm(j,i+nres)-scala2*
443      &          dc_norm(j,i))/vbld(i+1)-0.25d0*cosa/fac11*
444      &          dcostheta(j,2,i+1)-fac11*dcosalpha(j,2,i))*fac12+
445      &          (scala2-fac11*cosa)*(0.25d0*sina/fac10*
446      &          dcostheta(j,2,i+1)+fac10*cosa/sina*dcosalpha(j,2,i)
447      &          ))/fac13
448                 domega(j,2,i)=-1/sino*dcosomega(j,2,i)          
449                 dcosomega(j,3,i)=1/fac10*((1/vbld(i+nres)*(dc_norm(j,i)-
450      &          scala2*dc_norm(j,i+nres))-fac11*dcosalpha(j,3,i))*sina+
451      &          (scala2-fac11*cosa)*(cosa/sina*dcosalpha(j,3,i)))/fac14
452                 domega(j,3,i)=-1/sino*dcosomega(j,3,i)                          
453             enddo           
454           endif
455         endif    
456        enddo                                          
457 #endif
458 #if defined(MPI) && defined(PARINTDER)
459       if (nfgtasks.gt.1) then
460 #ifdef DEBUG
461        write (iout,*) "Gather dtheta"
462 cd      call flush(iout)
463 c      write (iout,*) "dtheta before gather"
464 c      do i=1,nres
465 c        write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),k=1,3),j=1,2)
466 c      enddo
467 #endif
468       call MPI_Gatherv(dtheta(1,1,ithet_start),ithet_count(fg_rank),
469      &  MPI_THET,dtheta(1,1,1),ithet_count(0),ithet_displ(0),MPI_THET,
470      &  king,FG_COMM,IERROR)
471 #ifdef DEBUG
472 cd      write (iout,*) "Gather dphi"
473 cd      call flush(iout)
474       write (iout,*) "dphi before gather"
475       do i=1,nres
476         write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),k=1,3),j=1,3)
477       enddo
478 #endif
479       call MPI_Gatherv(dphi(1,1,iphi1_start),iphi1_count(fg_rank),
480      &  MPI_GAM,dphi(1,1,1),iphi1_count(0),iphi1_displ(0),MPI_GAM,
481      &  king,FG_COMM,IERROR)
482 cd      write (iout,*) "Gather dalpha"
483 cd      call flush(iout)
484 #ifdef CRYST_SC
485       call MPI_Gatherv(dalpha(1,1,ibond_start),ibond_count(fg_rank),
486      &  MPI_GAM,dalpha(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
487      &  king,FG_COMM,IERROR)
488 cd      write (iout,*) "Gather domega"
489 cd      call flush(iout)
490       call MPI_Gatherv(domega(1,1,ibond_start),ibond_count(fg_rank),
491      &  MPI_GAM,domega(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
492      &  king,FG_COMM,IERROR)
493 #endif
494       endif
495 #endif
496 #ifdef DEBUG
497       write (iout,*) "dtheta after gather"
498       do i=1,nres
499         write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),j=1,3),j=1,2)
500       enddo
501       write (iout,*) "dphi after gather"
502       do i=1,nres
503         write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),j=1,3),k=1,3)
504       enddo
505 #endif
506       return
507       end
508        
509       subroutine checkintcartgrad
510       implicit real*8 (a-h,o-z)
511       include 'DIMENSIONS'
512 #ifdef MPI
513       include 'mpif.h'
514 #endif
515       include 'COMMON.CHAIN' 
516       include 'COMMON.VAR'
517       include 'COMMON.GEO'
518       include 'COMMON.INTERACT'
519       include 'COMMON.DERIV'
520       include 'COMMON.IOUNITS'
521       include 'COMMON.SETUP'
522       double precision dthetanum(3,2,maxres),dphinum(3,3,maxres)
523      & ,dalphanum(3,3,maxres), domeganum(3,3,maxres)
524       double precision theta_s(maxres),phi_s(maxres),alph_s(maxres),
525      & omeg_s(maxres),dc_norm_s(3)
526       double precision aincr /1.0d-5/
527       
528       do i=1,nres
529         phi_s(i)=phi(i)
530         theta_s(i)=theta(i)     
531         alph_s(i)=alph(i)
532         omeg_s(i)=omeg(i)
533       enddo
534 c Check theta gradient
535       write (iout,*) 
536      & "Analytical (upper) and numerical (lower) gradient of theta"
537       write (iout,*) 
538       do i=3,nres
539         do j=1,3
540           dcji=dc(j,i-2)
541           dc(j,i-2)=dcji+aincr
542           call chainbuild_cart
543           call int_from_cart1(.false.)
544           dthetanum(j,1,i)=(theta(i)-theta_s(i))/aincr 
545           dc(j,i-2)=dcji
546           dcji=dc(j,i-1)
547           dc(j,i-1)=dc(j,i-1)+aincr
548           call chainbuild_cart    
549           dthetanum(j,2,i)=(theta(i)-theta_s(i))/aincr
550           dc(j,i-1)=dcji
551         enddo 
552         write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dtheta(j,1,i),j=1,3),
553      &    (dtheta(j,2,i),j=1,3)
554         write (iout,'(5x,3f10.5,5x,3f10.5)') (dthetanum(j,1,i),j=1,3),
555      &    (dthetanum(j,2,i),j=1,3)
556         write (iout,'(5x,3f10.5,5x,3f10.5)') 
557      &    (dthetanum(j,1,i)/dtheta(j,1,i),j=1,3),
558      &    (dthetanum(j,2,i)/dtheta(j,2,i),j=1,3)
559         write (iout,*)
560       enddo
561 c Check gamma gradient
562       write (iout,*) 
563      & "Analytical (upper) and numerical (lower) gradient of gamma"
564       do i=4,nres
565         do j=1,3
566           dcji=dc(j,i-3)
567           dc(j,i-3)=dcji+aincr
568           call chainbuild_cart
569           dphinum(j,1,i)=(phi(i)-phi_s(i))/aincr  
570           dc(j,i-3)=dcji
571           dcji=dc(j,i-2)
572           dc(j,i-2)=dcji+aincr
573           call chainbuild_cart
574           dphinum(j,2,i)=(phi(i)-phi_s(i))/aincr 
575           dc(j,i-2)=dcji
576           dcji=dc(j,i-1)
577           dc(j,i-1)=dc(j,i-1)+aincr
578           call chainbuild_cart
579           dphinum(j,3,i)=(phi(i)-phi_s(i))/aincr
580           dc(j,i-1)=dcji
581         enddo 
582         write (iout,'(i5,3(3f10.5,5x))') i,(dphi(j,1,i),j=1,3),
583      &    (dphi(j,2,i),j=1,3),(dphi(j,3,i),j=1,3)
584         write (iout,'(5x,3(3f10.5,5x))') (dphinum(j,1,i),j=1,3),
585      &    (dphinum(j,2,i),j=1,3),(dphinum(j,3,i),j=1,3)
586         write (iout,'(5x,3(3f10.5,5x))') 
587      &    (dphinum(j,1,i)/dphi(j,1,i),j=1,3),
588      &    (dphinum(j,2,i)/dphi(j,2,i),j=1,3),
589      &    (dphinum(j,3,i)/dphi(j,3,i),j=1,3)
590         write (iout,*)
591       enddo
592 c Check alpha gradient
593       write (iout,*) 
594      & "Analytical (upper) and numerical (lower) gradient of alpha"
595       do i=2,nres-1
596        if(itype(i).ne.10) then
597             do j=1,3
598               dcji=dc(j,i-1)
599               dc(j,i-1)=dcji+aincr
600               call chainbuild_cart
601               dalphanum(j,1,i)=(alph(i)-alph_s(i))
602      &        /aincr  
603               dc(j,i-1)=dcji
604               dcji=dc(j,i)
605               dc(j,i)=dcji+aincr
606               call chainbuild_cart
607               dalphanum(j,2,i)=(alph(i)-alph_s(i))
608      &        /aincr 
609               dc(j,i)=dcji
610               dcji=dc(j,i+nres)
611               dc(j,i+nres)=dc(j,i+nres)+aincr
612               call chainbuild_cart
613               dalphanum(j,3,i)=(alph(i)-alph_s(i))
614      &        /aincr
615              dc(j,i+nres)=dcji
616             enddo
617           endif      
618         write (iout,'(i5,3(3f10.5,5x))') i,(dalpha(j,1,i),j=1,3),
619      &    (dalpha(j,2,i),j=1,3),(dalpha(j,3,i),j=1,3)
620         write (iout,'(5x,3(3f10.5,5x))') (dalphanum(j,1,i),j=1,3),
621      &    (dalphanum(j,2,i),j=1,3),(dalphanum(j,3,i),j=1,3)
622         write (iout,'(5x,3(3f10.5,5x))') 
623      &    (dalphanum(j,1,i)/dalpha(j,1,i),j=1,3),
624      &    (dalphanum(j,2,i)/dalpha(j,2,i),j=1,3),
625      &    (dalphanum(j,3,i)/dalpha(j,3,i),j=1,3)
626         write (iout,*)
627       enddo
628 c     Check omega gradient
629       write (iout,*) 
630      & "Analytical (upper) and numerical (lower) gradient of omega"
631       do i=2,nres-1
632        if(itype(i).ne.10) then
633             do j=1,3
634               dcji=dc(j,i-1)
635               dc(j,i-1)=dcji+aincr
636               call chainbuild_cart
637               domeganum(j,1,i)=(omeg(i)-omeg_s(i))
638      &        /aincr  
639               dc(j,i-1)=dcji
640               dcji=dc(j,i)
641               dc(j,i)=dcji+aincr
642               call chainbuild_cart
643               domeganum(j,2,i)=(omeg(i)-omeg_s(i))
644      &        /aincr 
645               dc(j,i)=dcji
646               dcji=dc(j,i+nres)
647               dc(j,i+nres)=dc(j,i+nres)+aincr
648               call chainbuild_cart
649               domeganum(j,3,i)=(omeg(i)-omeg_s(i))
650      &        /aincr
651              dc(j,i+nres)=dcji
652             enddo
653           endif      
654         write (iout,'(i5,3(3f10.5,5x))') i,(domega(j,1,i),j=1,3),
655      &    (domega(j,2,i),j=1,3),(domega(j,3,i),j=1,3)
656         write (iout,'(5x,3(3f10.5,5x))') (domeganum(j,1,i),j=1,3),
657      &    (domeganum(j,2,i),j=1,3),(domeganum(j,3,i),j=1,3)
658         write (iout,'(5x,3(3f10.5,5x))') 
659      &    (domeganum(j,1,i)/domega(j,1,i),j=1,3),
660      &    (domeganum(j,2,i)/domega(j,2,i),j=1,3),
661      &    (domeganum(j,3,i)/domega(j,3,i),j=1,3)
662         write (iout,*)
663       enddo
664       return
665       end
666
667       subroutine chainbuild_cart
668       implicit real*8 (a-h,o-z)
669       include 'DIMENSIONS'
670 #ifdef MPI
671       include 'mpif.h'
672 #endif
673       include 'COMMON.SETUP'
674       include 'COMMON.CHAIN' 
675       include 'COMMON.LOCAL'
676       include 'COMMON.TIME1'
677       include 'COMMON.IOUNITS'
678       
679 #ifdef MPI
680       if (nfgtasks.gt.1) then
681 c        write (iout,*) "BCAST in chainbuild_cart"
682 c        call flush(iout)
683 c Broadcast the order to build the chain and compute internal coordinates
684 c to the slaves. The slaves receive the order in ERGASTULUM.
685         time00=MPI_Wtime()
686 c      write (iout,*) "CHAINBUILD_CART: DC before BCAST"
687 c      do i=0,nres
688 c        write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
689 c     &   (dc(j,i+nres),j=1,3)
690 c      enddo 
691         if (fg_rank.eq.0) 
692      &    call MPI_Bcast(7,1,MPI_INTEGER,king,FG_COMM,IERROR)
693         time_bcast7=time_bcast7+MPI_Wtime()-time00
694         time01=MPI_Wtime()
695         call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,
696      &    king,FG_COMM,IERR)
697 c      write (iout,*) "CHAINBUILD_CART: DC after BCAST"
698 c      do i=0,nres
699 c        write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
700 c     &   (dc(j,i+nres),j=1,3)
701 c      enddo 
702 c        write (iout,*) "End BCAST in chainbuild_cart"
703 c        call flush(iout)
704         time_bcast=time_bcast+MPI_Wtime()-time00
705         time_bcastc=time_bcastc+MPI_Wtime()-time01
706       endif
707 #endif
708       do j=1,3
709         c(j,1)=dc(j,0)
710       enddo
711       do i=2,nres
712         do j=1,3
713           c(j,i)=c(j,i-1)+dc(j,i-1)
714         enddo
715       enddo 
716       do i=1,nres
717         do j=1,3
718           c(j,i+nres)=c(j,i)+dc(j,i+nres)
719         enddo
720       enddo
721 c      write (iout,*) "CHAINBUILD_CART"
722 c      call cartprint
723       call int_from_cart1(.false.)
724       return
725       end