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