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