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