ctest wham & cluster src-M
[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).ge.-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       do i=1,nres-1
153        do j=1,3
154         dc_norm2(j,i+nres)=-dc_norm(j,i+nres)
155        enddo
156       enddo
157 Calculate derivative of Tauangle
158 #ifdef PARINTDER
159       do i=itau_start,itau_end
160 #else
161       do i=3,nres
162 #endif
163        if ((itype(i-2).eq.21).or.(itype(i-2).eq.10)) cycle
164 cc dtauangle(j,intertyp,dervityp,residue number)
165 cc INTERTYP=1 SC...Ca...Ca..Ca
166 c the conventional case
167         sint=dsin(theta(i))
168         sint1=dsin(omicron(2,i-1))
169         sing=dsin(tauangle(1,i))
170         cost=dcos(theta(i))
171         cost1=dcos(omicron(2,i-1))
172         cosg=dcos(tauangle(1,i))
173 C        do j=1,3
174 C        dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
175 cc       write(iout,*) dc_norm2(j,i-2+nres),"dcnorm"
176 C        enddo
177         scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1))
178         fac0=1.0d0/(sint1*sint)
179         fac1=cost*fac0
180         fac2=cost1*fac0
181         fac3=cosg*cost1/(sint1*sint1)
182         fac4=cosg*cost/(sint*sint)
183 cc         write(iout,*) "faki",fac0,fac1,fac2,fac3,fac4
184 c    Obtaining the gamma derivatives from sine derivative                                
185        if (tauangle(1,i).gt.-pi4.and.tauangle(1,i).le.pi4.or.
186      &     tauangle(1,i).gt.pi34.and.tauangle(1,i).le.pi.or.
187      &     tauangle(1,i).gt.-pi.and.tauangle(1,i).le.-pi34) then
188          call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
189          call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1),vp2)
190          call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
191         do j=1,3
192             ctgt=cost/sint
193             ctgt1=cost1/sint1
194             cosg_inv=1.0d0/cosg
195             dsintau(j,1,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
196      &-(fac0*vp1(j)+sing*(dc_norm2(j,i-2+nres)))
197      & *vbld_inv(i-2+nres)
198             dtauangle(j,1,1,i)=cosg_inv*dsintau(j,1,1,i)
199             dsintau(j,1,2,i)=
200      &        -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*dtheta(j,1,i))
201      &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
202 c            write(iout,*) "dsintau", dsintau(j,1,1,i),dsintau(j,1,2,i)
203             dtauangle(j,1,2,i)=cosg_inv*dsintau(j,1,2,i)
204 c Bug fixed 3/24/05 (AL)
205             dsintau(j,1,3,i)=-sing*ctgt*dtheta(j,2,i)
206      &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
207 c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
208             dtauangle(j,1,3,i)=cosg_inv*dsintau(j,1,3,i)
209          enddo                                          
210 c   Obtaining the gamma derivatives from cosine derivative
211         else
212            do j=1,3
213            dcostau(j,1,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
214      &     dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1)-scalp*
215      &     (dc_norm2(j,i-2+nres)))/vbld(i-2+nres)
216            dtauangle(j,1,1,i)=-1/sing*dcostau(j,1,1,i)
217            dcostau(j,1,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
218      &     dcostheta(j,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
219      &     dcostheta(j,1,i)
220            dtauangle(j,1,2,i)=-1/sing*dcostau(j,1,2,i)
221            dcostau(j,1,3,i)=fac2*dcostheta(j,2,i)+fac4*
222      &     dcostheta(j,2,i)-fac0*(-dc_norm(j,i-2+nres)-scalp*
223      &     dc_norm(j,i-1))/vbld(i)
224            dtauangle(j,1,3,i)=-1/sing*dcostau(j,1,3,i)
225 c         write (iout,*) "else",i
226          enddo
227         endif
228 c        do k=1,3                 
229 c        write(iout,*) "tu",1,i,k,(dtauangle(j,1,k,i),j=1,3)        
230 c        write(iout,*) "tu",2,i,k,(dtauangle(j,2,k,i),j=1,3)        
231 c        write(iout,*) "tu",3,i,k,(dtauangle(j,3,k,i),j=1,3)        
232 c        enddo                
233       enddo
234 CC Second case Ca...Ca...Ca...SC
235 #ifdef PARINTDER
236       do i=itau_start,itau_end
237 #else
238       do i=4,nres
239 #endif
240        if ((itype(i-1).eq.21).or.(itype(i-1).eq.10)) cycle
241 c the conventional case
242         sint=dsin(omicron(1,i))
243         sint1=dsin(theta(i-1))
244         sing=dsin(tauangle(2,i))
245         cost=dcos(omicron(1,i))
246         cost1=dcos(theta(i-1))
247         cosg=dcos(tauangle(2,i))
248 c        do j=1,3
249 c        dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
250 c        enddo
251         scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1+nres))
252         fac0=1.0d0/(sint1*sint)
253         fac1=cost*fac0
254         fac2=cost1*fac0
255         fac3=cosg*cost1/(sint1*sint1)
256         fac4=cosg*cost/(sint*sint)
257 c    Obtaining the gamma derivatives from sine derivative                                
258        if (tauangle(2,i).gt.-pi4.and.tauangle(2,i).le.pi4.or.
259      &     tauangle(2,i).gt.pi34.and.tauangle(2,i).le.pi.or.
260      &     tauangle(2,i).gt.-pi.and.tauangle(2,i).le.-pi34) then
261          call vecpr(dc_norm2(1,i-1+nres),dc_norm(1,i-2),vp1)
262          call vecpr(dc_norm(1,i-3),dc_norm(1,i-1+nres),vp2)
263          call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
264         do j=1,3
265             ctgt=cost/sint
266             ctgt1=cost1/sint1
267             cosg_inv=1.0d0/cosg
268             dsintau(j,2,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
269      &        +(fac0*vp1(j)-sing*dc_norm(j,i-3))*vbld_inv(i-2)
270 c       write(iout,*) i,j,dsintau(j,2,1,i),sing*ctgt1*dtheta(j,1,i-1),
271 c     &fac0*vp1(j),sing*dc_norm(j,i-3),vbld_inv(i-2),"dsintau(2,1)"
272             dtauangle(j,2,1,i)=cosg_inv*dsintau(j,2,1,i)
273             dsintau(j,2,2,i)=
274      &        -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*domicron(j,1,1,i))
275      &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
276 c            write(iout,*) "sprawdzenie",i,j,sing*ctgt1*dtheta(j,2,i-1),
277 c     & sing*ctgt*domicron(j,1,2,i),
278 c     & (fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
279             dtauangle(j,2,2,i)=cosg_inv*dsintau(j,2,2,i)
280 c Bug fixed 3/24/05 (AL)
281             dsintau(j,2,3,i)=-sing*ctgt*domicron(j,1,2,i)
282      &       +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))*vbld_inv(i-1+nres)
283 c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
284             dtauangle(j,2,3,i)=cosg_inv*dsintau(j,2,3,i)
285          enddo                                          
286 c   Obtaining the gamma derivatives from cosine derivative
287         else
288            do j=1,3
289            dcostau(j,2,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
290      &     dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
291      &     dc_norm(j,i-3))/vbld(i-2)
292            dtauangle(j,2,1,i)=-1/sing*dcostau(j,2,1,i)
293            dcostau(j,2,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
294      &     dcosomicron(j,1,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
295      &     dcosomicron(j,1,1,i)
296            dtauangle(j,2,2,i)=-1/sing*dcostau(j,2,2,i)
297            dcostau(j,2,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
298      &     dcosomicron(j,1,2,i)-fac0*(dc_norm(j,i-3)-scalp*
299      &     dc_norm(j,i-1+nres))/vbld(i-1+nres)
300            dtauangle(j,2,3,i)=-1/sing*dcostau(j,2,3,i)
301 c        write(iout,*) i,j,"else", dtauangle(j,2,3,i) 
302          enddo
303         endif                                                                                            
304       enddo
305
306
307 CCC third case SC...Ca...Ca...SC
308 #ifdef PARINTDER
309
310       do i=itau_start,itau_end
311 #else
312       do i=3,nres
313 #endif
314 c the conventional case
315       if ((itype(i-1).eq.21).or.(itype(i-1).eq.10).or.
316      &(itype(i-2).eq.21).or.(itype(i-2).eq.10)) cycle
317         sint=dsin(omicron(1,i))
318         sint1=dsin(omicron(2,i-1))
319         sing=dsin(tauangle(3,i))
320         cost=dcos(omicron(1,i))
321         cost1=dcos(omicron(2,i-1))
322         cosg=dcos(tauangle(3,i))
323 C        do j=1,3
324 C        dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
325 c        dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
326 C        enddo
327         scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres))
328         fac0=1.0d0/(sint1*sint)
329         fac1=cost*fac0
330         fac2=cost1*fac0
331         fac3=cosg*cost1/(sint1*sint1)
332         fac4=cosg*cost/(sint*sint)
333 c    Obtaining the gamma derivatives from sine derivative                                
334        if (tauangle(3,i).gt.-pi4.and.tauangle(3,i).le.pi4.or.
335      &     tauangle(3,i).gt.pi34.and.tauangle(3,i).le.pi.or.
336      &     tauangle(3,i).gt.-pi.and.tauangle(3,i).le.-pi34) then
337          call vecpr(dc_norm(1,i-1+nres),dc_norm(1,i-2),vp1)
338          call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres),vp2)
339          call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
340         do j=1,3
341             ctgt=cost/sint
342             ctgt1=cost1/sint1
343             cosg_inv=1.0d0/cosg
344             dsintau(j,3,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
345      &        -(fac0*vp1(j)-sing*dc_norm(j,i-2+nres))
346      &        *vbld_inv(i-2+nres)
347             dtauangle(j,3,1,i)=cosg_inv*dsintau(j,3,1,i)
348             dsintau(j,3,2,i)=
349      &        -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*domicron(j,1,1,i))
350      &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
351             dtauangle(j,3,2,i)=cosg_inv*dsintau(j,3,2,i)
352 c Bug fixed 3/24/05 (AL)
353             dsintau(j,3,3,i)=-sing*ctgt*domicron(j,1,2,i)
354      &        +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))
355      &        *vbld_inv(i-1+nres)
356 c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
357             dtauangle(j,3,3,i)=cosg_inv*dsintau(j,3,3,i)
358          enddo                                          
359 c   Obtaining the gamma derivatives from cosine derivative
360         else
361            do j=1,3
362            dcostau(j,3,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
363      &     dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
364      &     dc_norm2(j,i-2+nres))/vbld(i-2+nres)
365            dtauangle(j,3,1,i)=-1/sing*dcostau(j,3,1,i)
366            dcostau(j,3,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
367      &     dcosomicron(j,1,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
368      &     dcosomicron(j,1,1,i)
369            dtauangle(j,3,2,i)=-1/sing*dcostau(j,3,2,i)
370            dcostau(j,3,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
371      &     dcosomicron(j,1,2,i)-fac0*(dc_norm2(j,i-2+nres)-scalp*
372      &     dc_norm(j,i-1+nres))/vbld(i-1+nres)
373            dtauangle(j,3,3,i)=-1/sing*dcostau(j,3,3,i)
374 c          write(iout,*) "else",i 
375          enddo
376         endif                                                                                            
377       enddo
378 #ifdef CRYST_SC
379 c   Derivatives of side-chain angles alpha and omega
380 #if defined(MPI) && defined(PARINTDER)
381         do i=ibond_start,ibond_end
382 #else
383         do i=2,nres-1           
384 #endif
385           if(itype(i).ne.10) then         
386              fac5=1.0d0/dsqrt(2*(1+dcos(theta(i+1))))
387              fac6=fac5/vbld(i)
388              fac7=fac5*fac5
389              fac8=fac5/vbld(i+1)     
390              fac9=fac5/vbld(i+nres)                  
391              scala1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
392              scala2=scalar(dc_norm(1,i),dc_norm(1,i+nres))
393              cosa=dsqrt(0.5d0/(1.0d0+dcos(theta(i+1))))*(
394      &       scalar(dC_norm(1,i),dC_norm(1,i+nres))
395      &       -scalar(dC_norm(1,i-1),dC_norm(1,i+nres)))
396              sina=sqrt(1-cosa*cosa)
397              sino=dsin(omeg(i))                                                                                              
398              do j=1,3     
399                 dcosalpha(j,1,i)=fac6*(scala1*dc_norm(j,i-1)-
400      &          dc_norm(j,i+nres))-cosa*fac7*dcostheta(j,1,i+1)
401                 dalpha(j,1,i)=-1/sina*dcosalpha(j,1,i)
402                 dcosalpha(j,2,i)=fac8*(dc_norm(j,i+nres)-
403      &          scala2*dc_norm(j,i))-cosa*fac7*dcostheta(j,2,i+1)
404                 dalpha(j,2,i)=-1/sina*dcosalpha(j,2,i)
405                 dcosalpha(j,3,i)=(fac9*(dc_norm(j,i)-
406      &          dc_norm(j,i-1))-(cosa*dc_norm(j,i+nres))/
407      &          vbld(i+nres))
408                 dalpha(j,3,i)=-1/sina*dcosalpha(j,3,i)
409             enddo
410 c obtaining the derivatives of omega from sines     
411             if(omeg(i).gt.-pi4.and.omeg(i).le.pi4.or.
412      &         omeg(i).gt.pi34.and.omeg(i).le.pi.or.
413      &         omeg(i).gt.-pi.and.omeg(i).le.-pi34) then
414                fac15=dcos(theta(i+1))/(dsin(theta(i+1))*
415      &         dsin(theta(i+1)))
416                fac16=dcos(alph(i))/(dsin(alph(i))*dsin(alph(i)))
417                fac17=1.0d0/(dsin(theta(i+1))*dsin(alph(i)))             
418                call vecpr(dc_norm(1,i+nres),dc_norm(1,i),vo1)
419                call vecpr(dc_norm(1,i+nres),dc_norm(1,i-1),vo2)
420                call vecpr(dc_norm(1,i),dc_norm(1,i-1),vo3)
421                coso_inv=1.0d0/dcos(omeg(i))                            
422                do j=1,3
423                  dsinomega(j,1,i)=sino*(fac15*dcostheta(j,1,i+1)
424      &           +fac16*dcosalpha(j,1,i))-fac17/vbld(i)*vo1(j)-(
425      &           sino*dc_norm(j,i-1))/vbld(i)
426                  domega(j,1,i)=coso_inv*dsinomega(j,1,i)
427                  dsinomega(j,2,i)=sino*(fac15*dcostheta(j,2,i+1)
428      &           +fac16*dcosalpha(j,2,i))+fac17/vbld(i+1)*vo2(j)
429      &           -sino*dc_norm(j,i)/vbld(i+1)
430                  domega(j,2,i)=coso_inv*dsinomega(j,2,i)                                                       
431                  dsinomega(j,3,i)=sino*fac16*dcosalpha(j,3,i)-
432      &           fac17/vbld(i+nres)*vo3(j)-sino*dc_norm(j,i+nres)/
433      &           vbld(i+nres)
434                  domega(j,3,i)=coso_inv*dsinomega(j,3,i)
435               enddo                              
436            else
437 c   obtaining the derivatives of omega from cosines
438              fac10=sqrt(0.5d0*(1-dcos(theta(i+1))))
439              fac11=sqrt(0.5d0*(1+dcos(theta(i+1))))
440              fac12=fac10*sina
441              fac13=fac12*fac12
442              fac14=sina*sina
443              do j=1,3                                    
444                 dcosomega(j,1,i)=(-(0.25d0*cosa/fac11*
445      &          dcostheta(j,1,i+1)+fac11*dcosalpha(j,1,i))*fac12+
446      &          (0.25d0/fac10*sina*dcostheta(j,1,i+1)+cosa/sina*
447      &          fac10*dcosalpha(j,1,i))*(scala2-fac11*cosa))/fac13
448                 domega(j,1,i)=-1/sino*dcosomega(j,1,i)
449                 dcosomega(j,2,i)=(((dc_norm(j,i+nres)-scala2*
450      &          dc_norm(j,i))/vbld(i+1)-0.25d0*cosa/fac11*
451      &          dcostheta(j,2,i+1)-fac11*dcosalpha(j,2,i))*fac12+
452      &          (scala2-fac11*cosa)*(0.25d0*sina/fac10*
453      &          dcostheta(j,2,i+1)+fac10*cosa/sina*dcosalpha(j,2,i)
454      &          ))/fac13
455                 domega(j,2,i)=-1/sino*dcosomega(j,2,i)          
456                 dcosomega(j,3,i)=1/fac10*((1/vbld(i+nres)*(dc_norm(j,i)-
457      &          scala2*dc_norm(j,i+nres))-fac11*dcosalpha(j,3,i))*sina+
458      &          (scala2-fac11*cosa)*(cosa/sina*dcosalpha(j,3,i)))/fac14
459                 domega(j,3,i)=-1/sino*dcosomega(j,3,i)                          
460             enddo           
461           endif
462         endif    
463        enddo                                          
464 #endif
465 #if defined(MPI) && defined(PARINTDER)
466       if (nfgtasks.gt.1) then
467 #ifdef DEBUG
468        write (iout,*) "Gather dtheta"
469 cd      call flush(iout)
470 c      write (iout,*) "dtheta before gather"
471 c      do i=1,nres
472 c        write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),k=1,3),j=1,2)
473 c      enddo
474 #endif
475       call MPI_Gatherv(dtheta(1,1,ithet_start),ithet_count(fg_rank),
476      &  MPI_THET,dtheta(1,1,1),ithet_count(0),ithet_displ(0),MPI_THET,
477      &  king,FG_COMM,IERROR)
478 #ifdef DEBUG
479 cd      write (iout,*) "Gather dphi"
480 cd      call flush(iout)
481       write (iout,*) "dphi before gather"
482       do i=1,nres
483         write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),k=1,3),j=1,3)
484       enddo
485 #endif
486       call MPI_Gatherv(dphi(1,1,iphi1_start),iphi1_count(fg_rank),
487      &  MPI_GAM,dphi(1,1,1),iphi1_count(0),iphi1_displ(0),MPI_GAM,
488      &  king,FG_COMM,IERROR)
489 cd      write (iout,*) "Gather dalpha"
490 cd      call flush(iout)
491 #ifdef CRYST_SC
492       call MPI_Gatherv(dalpha(1,1,ibond_start),ibond_count(fg_rank),
493      &  MPI_GAM,dalpha(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
494      &  king,FG_COMM,IERROR)
495 cd      write (iout,*) "Gather domega"
496 cd      call flush(iout)
497       call MPI_Gatherv(domega(1,1,ibond_start),ibond_count(fg_rank),
498      &  MPI_GAM,domega(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
499      &  king,FG_COMM,IERROR)
500 #endif
501       endif
502 #endif
503 #ifdef DEBUG
504       write (iout,*) "dtheta after gather"
505       do i=1,nres
506         write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),j=1,3),j=1,2)
507       enddo
508       write (iout,*) "dphi after gather"
509       do i=1,nres
510         write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),j=1,3),k=1,3)
511       enddo
512 #endif
513       return
514       end
515        
516       subroutine checkintcartgrad
517       implicit real*8 (a-h,o-z)
518       include 'DIMENSIONS'
519 #ifdef MPI
520       include 'mpif.h'
521 #endif
522       include 'COMMON.CHAIN' 
523       include 'COMMON.VAR'
524       include 'COMMON.GEO'
525       include 'COMMON.INTERACT'
526       include 'COMMON.DERIV'
527       include 'COMMON.IOUNITS'
528       include 'COMMON.SETUP'
529       double precision dthetanum(3,2,maxres),dphinum(3,3,maxres)
530      & ,dalphanum(3,3,maxres), domeganum(3,3,maxres)
531       double precision theta_s(maxres),phi_s(maxres),alph_s(maxres),
532      & omeg_s(maxres),dc_norm_s(3)
533       double precision aincr /1.0d-5/
534       
535       do i=1,nres
536         phi_s(i)=phi(i)
537         theta_s(i)=theta(i)     
538         alph_s(i)=alph(i)
539         omeg_s(i)=omeg(i)
540       enddo
541 c Check theta gradient
542       write (iout,*) 
543      & "Analytical (upper) and numerical (lower) gradient of theta"
544       write (iout,*) 
545       do i=3,nres
546         do j=1,3
547           dcji=dc(j,i-2)
548           dc(j,i-2)=dcji+aincr
549           call chainbuild_cart
550           call int_from_cart1(.false.)
551           dthetanum(j,1,i)=(theta(i)-theta_s(i))/aincr 
552           dc(j,i-2)=dcji
553           dcji=dc(j,i-1)
554           dc(j,i-1)=dc(j,i-1)+aincr
555           call chainbuild_cart    
556           dthetanum(j,2,i)=(theta(i)-theta_s(i))/aincr
557           dc(j,i-1)=dcji
558         enddo 
559         write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dtheta(j,1,i),j=1,3),
560      &    (dtheta(j,2,i),j=1,3)
561         write (iout,'(5x,3f10.5,5x,3f10.5)') (dthetanum(j,1,i),j=1,3),
562      &    (dthetanum(j,2,i),j=1,3)
563         write (iout,'(5x,3f10.5,5x,3f10.5)') 
564      &    (dthetanum(j,1,i)/dtheta(j,1,i),j=1,3),
565      &    (dthetanum(j,2,i)/dtheta(j,2,i),j=1,3)
566         write (iout,*)
567       enddo
568 c Check gamma gradient
569       write (iout,*) 
570      & "Analytical (upper) and numerical (lower) gradient of gamma"
571       do i=4,nres
572         do j=1,3
573           dcji=dc(j,i-3)
574           dc(j,i-3)=dcji+aincr
575           call chainbuild_cart
576           dphinum(j,1,i)=(phi(i)-phi_s(i))/aincr  
577           dc(j,i-3)=dcji
578           dcji=dc(j,i-2)
579           dc(j,i-2)=dcji+aincr
580           call chainbuild_cart
581           dphinum(j,2,i)=(phi(i)-phi_s(i))/aincr 
582           dc(j,i-2)=dcji
583           dcji=dc(j,i-1)
584           dc(j,i-1)=dc(j,i-1)+aincr
585           call chainbuild_cart
586           dphinum(j,3,i)=(phi(i)-phi_s(i))/aincr
587           dc(j,i-1)=dcji
588         enddo 
589         write (iout,'(i5,3(3f10.5,5x))') i,(dphi(j,1,i),j=1,3),
590      &    (dphi(j,2,i),j=1,3),(dphi(j,3,i),j=1,3)
591         write (iout,'(5x,3(3f10.5,5x))') (dphinum(j,1,i),j=1,3),
592      &    (dphinum(j,2,i),j=1,3),(dphinum(j,3,i),j=1,3)
593         write (iout,'(5x,3(3f10.5,5x))') 
594      &    (dphinum(j,1,i)/dphi(j,1,i),j=1,3),
595      &    (dphinum(j,2,i)/dphi(j,2,i),j=1,3),
596      &    (dphinum(j,3,i)/dphi(j,3,i),j=1,3)
597         write (iout,*)
598       enddo
599 c Check alpha gradient
600       write (iout,*) 
601      & "Analytical (upper) and numerical (lower) gradient of alpha"
602       do i=2,nres-1
603        if(itype(i).ne.10) then
604             do j=1,3
605               dcji=dc(j,i-1)
606               dc(j,i-1)=dcji+aincr
607               call chainbuild_cart
608               dalphanum(j,1,i)=(alph(i)-alph_s(i))
609      &        /aincr  
610               dc(j,i-1)=dcji
611               dcji=dc(j,i)
612               dc(j,i)=dcji+aincr
613               call chainbuild_cart
614               dalphanum(j,2,i)=(alph(i)-alph_s(i))
615      &        /aincr 
616               dc(j,i)=dcji
617               dcji=dc(j,i+nres)
618               dc(j,i+nres)=dc(j,i+nres)+aincr
619               call chainbuild_cart
620               dalphanum(j,3,i)=(alph(i)-alph_s(i))
621      &        /aincr
622              dc(j,i+nres)=dcji
623             enddo
624           endif      
625         write (iout,'(i5,3(3f10.5,5x))') i,(dalpha(j,1,i),j=1,3),
626      &    (dalpha(j,2,i),j=1,3),(dalpha(j,3,i),j=1,3)
627         write (iout,'(5x,3(3f10.5,5x))') (dalphanum(j,1,i),j=1,3),
628      &    (dalphanum(j,2,i),j=1,3),(dalphanum(j,3,i),j=1,3)
629         write (iout,'(5x,3(3f10.5,5x))') 
630      &    (dalphanum(j,1,i)/dalpha(j,1,i),j=1,3),
631      &    (dalphanum(j,2,i)/dalpha(j,2,i),j=1,3),
632      &    (dalphanum(j,3,i)/dalpha(j,3,i),j=1,3)
633         write (iout,*)
634       enddo
635 c     Check omega gradient
636       write (iout,*) 
637      & "Analytical (upper) and numerical (lower) gradient of omega"
638       do i=2,nres-1
639        if(itype(i).ne.10) then
640             do j=1,3
641               dcji=dc(j,i-1)
642               dc(j,i-1)=dcji+aincr
643               call chainbuild_cart
644               domeganum(j,1,i)=(omeg(i)-omeg_s(i))
645      &        /aincr  
646               dc(j,i-1)=dcji
647               dcji=dc(j,i)
648               dc(j,i)=dcji+aincr
649               call chainbuild_cart
650               domeganum(j,2,i)=(omeg(i)-omeg_s(i))
651      &        /aincr 
652               dc(j,i)=dcji
653               dcji=dc(j,i+nres)
654               dc(j,i+nres)=dc(j,i+nres)+aincr
655               call chainbuild_cart
656               domeganum(j,3,i)=(omeg(i)-omeg_s(i))
657      &        /aincr
658              dc(j,i+nres)=dcji
659             enddo
660           endif      
661         write (iout,'(i5,3(3f10.5,5x))') i,(domega(j,1,i),j=1,3),
662      &    (domega(j,2,i),j=1,3),(domega(j,3,i),j=1,3)
663         write (iout,'(5x,3(3f10.5,5x))') (domeganum(j,1,i),j=1,3),
664      &    (domeganum(j,2,i),j=1,3),(domeganum(j,3,i),j=1,3)
665         write (iout,'(5x,3(3f10.5,5x))') 
666      &    (domeganum(j,1,i)/domega(j,1,i),j=1,3),
667      &    (domeganum(j,2,i)/domega(j,2,i),j=1,3),
668      &    (domeganum(j,3,i)/domega(j,3,i),j=1,3)
669         write (iout,*)
670       enddo
671       return
672       end
673
674       subroutine chainbuild_cart
675       implicit real*8 (a-h,o-z)
676       include 'DIMENSIONS'
677 #ifdef MPI
678       include 'mpif.h'
679 #endif
680       include 'COMMON.SETUP'
681       include 'COMMON.CHAIN' 
682       include 'COMMON.LOCAL'
683       include 'COMMON.TIME1'
684       include 'COMMON.IOUNITS'
685       
686 #ifdef MPI
687       if (nfgtasks.gt.1) then
688 c        write (iout,*) "BCAST in chainbuild_cart"
689 c        call flush(iout)
690 c Broadcast the order to build the chain and compute internal coordinates
691 c to the slaves. The slaves receive the order in ERGASTULUM.
692         time00=MPI_Wtime()
693 c      write (iout,*) "CHAINBUILD_CART: DC before BCAST"
694 c      do i=0,nres
695 c        write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
696 c     &   (dc(j,i+nres),j=1,3)
697 c      enddo 
698         if (fg_rank.eq.0) 
699      &    call MPI_Bcast(7,1,MPI_INTEGER,king,FG_COMM,IERROR)
700         time_bcast7=time_bcast7+MPI_Wtime()-time00
701         time01=MPI_Wtime()
702         call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,
703      &    king,FG_COMM,IERR)
704 c      write (iout,*) "CHAINBUILD_CART: DC after BCAST"
705 c      do i=0,nres
706 c        write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
707 c     &   (dc(j,i+nres),j=1,3)
708 c      enddo 
709 c        write (iout,*) "End BCAST in chainbuild_cart"
710 c        call flush(iout)
711         time_bcast=time_bcast+MPI_Wtime()-time00
712         time_bcastc=time_bcastc+MPI_Wtime()-time01
713       endif
714 #endif
715       do j=1,3
716         c(j,1)=dc(j,0)
717       enddo
718       do i=2,nres
719         do j=1,3
720           c(j,i)=c(j,i-1)+dc(j,i-1)
721         enddo
722       enddo 
723       do i=1,nres
724         do j=1,3
725           c(j,i+nres)=c(j,i)+dc(j,i+nres)
726         enddo
727       enddo
728 c      write (iout,*) "CHAINBUILD_CART"
729 c      call cartprint
730       call int_from_cart1(.false.)
731       return
732       end