wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / source / unres / src_MD-M / chainbuild.F
1       subroutine chainbuild
2
3 C Build the virtual polypeptide chain. Side-chain centroids are moveable.
4 C As of 2/17/95.
5 C
6       implicit real*8 (a-h,o-z)
7       include 'DIMENSIONS'
8       include 'COMMON.CHAIN'
9       include 'COMMON.LOCAL'
10       include 'COMMON.GEO'
11       include 'COMMON.VAR'
12       include 'COMMON.IOUNITS'
13       include 'COMMON.NAMES'
14       include 'COMMON.INTERACT'
15       double precision e1(3),e2(3),e3(3)
16       logical lprn,perbox,fail
17 C Set lprn=.true. for debugging
18       lprn = .false.
19       perbox=.false.
20       fail=.false.
21       call chainbuild_cart
22       return
23       end
24 C#ifdef DEBUG
25 C      if (perbox) then
26 C first three CA's and SC(2).
27 C
28       subroutine chainbuild_extconf
29
30 C Build the virtual polypeptide chain. Side-chain centroids are moveable.
31 C As of 2/17/95.
32 C
33       implicit real*8 (a-h,o-z)
34       include 'DIMENSIONS'
35       include 'COMMON.CHAIN'
36       include 'COMMON.LOCAL'
37       include 'COMMON.GEO'
38       include 'COMMON.VAR'
39       include 'COMMON.IOUNITS'
40       include 'COMMON.NAMES'
41       include 'COMMON.INTERACT'
42       double precision e1(3),e2(3),e3(3)
43       logical lprn,perbox,fail
44
45       call orig_frame
46 *
47 * Build the alpha-carbon chain.
48 *
49       do i=4,nres
50         call locate_next_res(i)
51       enddo     
52 C
53 C First and last SC must coincide with the corresponding CA.
54 C
55       do j=1,3
56         dc(j,nres+1)=0.0D0
57         dc_norm(j,nres+1)=0.0D0
58         dc(j,nres+nres)=0.0D0
59         dc_norm(j,nres+nres)=0.0D0
60         c(j,nres+1)=c(j,1)
61         c(j,nres+nres)=c(j,nres)
62       enddo
63 *
64 * Temporary diagnosis
65 *
66       if (lprn) then
67
68       call cartprint
69       write (iout,'(/a)') 'Recalculated internal coordinates'
70       do i=2,nres-1
71         do j=1,3
72           c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
73         enddo
74         be=0.0D0
75         if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i)
76         be1=rad2deg*beta(nres+i,i,maxres2,i+1)
77         alfai=0.0D0
78         if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
79         write (iout,1212) restyp(itype(i)),i,dist(i-1,i),
80      &  alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,maxres2),be1
81       enddo   
82  1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
83
84 C      endif
85       endif
86       return
87       end
88 C#endif
89 c-------------------------------------------------------------------------
90       subroutine orig_frame
91 C
92 C Define the origin and orientation of the coordinate system and locate 
93 C the first three atoms.
94 C
95       implicit real*8 (a-h,o-z)
96       include 'DIMENSIONS'
97       include 'COMMON.CHAIN'
98       include 'COMMON.LOCAL'
99       include 'COMMON.GEO'
100       include 'COMMON.VAR'
101       cost=dcos(theta(3))
102       sint=dsin(theta(3))
103       t(1,1,1)=-cost
104       t(1,2,1)=-sint 
105       t(1,3,1)= 0.0D0
106       t(2,1,1)=-sint
107       t(2,2,1)= cost
108       t(2,3,1)= 0.0D0
109       t(3,1,1)= 0.0D0
110       t(3,2,1)= 0.0D0
111       t(3,3,1)= 1.0D0
112       r(1,1,1)= 1.0D0
113       r(1,2,1)= 0.0D0
114       r(1,3,1)= 0.0D0
115       r(2,1,1)= 0.0D0
116       r(2,2,1)= 1.0D0
117       r(2,3,1)= 0.0D0
118       r(3,1,1)= 0.0D0
119       r(3,2,1)= 0.0D0
120       r(3,3,1)= 1.0D0
121       do i=1,3
122         do j=1,3
123           rt(i,j,1)=t(i,j,1)
124         enddo
125       enddo
126       do i=1,3
127         do j=1,3
128           prod(i,j,1)=0.0D0
129           prod(i,j,2)=t(i,j,1)
130         enddo
131         prod(i,i,1)=1.0D0
132       enddo   
133       c(1,1)=0.0D0
134       c(2,1)=0.0D0
135       c(3,1)=0.0D0
136       c(1,2)=vbld(2)
137       c(2,2)=0.0D0
138       c(3,2)=0.0D0
139       dc(1,0)=0.0d0
140       dc(2,0)=0.0D0
141       dc(3,0)=0.0D0
142       dc(1,1)=vbld(2)
143       dc(2,1)=0.0D0
144       dc(3,1)=0.0D0
145       dc_norm(1,0)=0.0D0
146       dc_norm(2,0)=0.0D0
147       dc_norm(3,0)=0.0D0
148       dc_norm(1,1)=1.0D0
149       dc_norm(2,1)=0.0D0
150       dc_norm(3,1)=0.0D0
151       do j=1,3
152         dc_norm(j,2)=prod(j,1,2)
153         dc(j,2)=vbld(3)*prod(j,1,2)
154         c(j,3)=c(j,2)+dc(j,2)
155       enddo
156       call locate_side_chain(2)
157       return
158       end
159 c-----------------------------------------------------------------------------
160       subroutine locate_next_res(i)
161 C
162 C Locate CA(i) and SC(i-1)
163 C
164       implicit real*8 (a-h,o-z)
165       include 'DIMENSIONS'
166       include 'COMMON.CHAIN'
167       include 'COMMON.LOCAL'
168       include 'COMMON.GEO'
169       include 'COMMON.VAR'
170       include 'COMMON.IOUNITS'
171       include 'COMMON.NAMES'
172       include 'COMMON.INTERACT'
173 C
174 C Define the rotation matrices corresponding to CA(i)
175 C
176 #ifdef OSF
177       theti=theta(i)
178       if (theti.ne.theti) theti=100.0     
179       phii=phi(i)
180       if (phii.ne.phii) phii=180.0     
181 #else
182       theti=theta(i)      
183       phii=phi(i)
184 #endif
185       cost=dcos(theti)
186       sint=dsin(theti)
187       cosphi=dcos(phii)
188       sinphi=dsin(phii)
189 * Define the matrices of the rotation about the virtual-bond valence angles
190 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
191 * program), R(i,j,k), and, the cumulative matrices of rotation RT
192       t(1,1,i-2)=-cost
193       t(1,2,i-2)=-sint 
194       t(1,3,i-2)= 0.0D0
195       t(2,1,i-2)=-sint
196       t(2,2,i-2)= cost
197       t(2,3,i-2)= 0.0D0
198       t(3,1,i-2)= 0.0D0
199       t(3,2,i-2)= 0.0D0
200       t(3,3,i-2)= 1.0D0
201       r(1,1,i-2)= 1.0D0
202       r(1,2,i-2)= 0.0D0
203       r(1,3,i-2)= 0.0D0
204       r(2,1,i-2)= 0.0D0
205       r(2,2,i-2)=-cosphi
206       r(2,3,i-2)= sinphi
207       r(3,1,i-2)= 0.0D0
208       r(3,2,i-2)= sinphi
209       r(3,3,i-2)= cosphi
210       rt(1,1,i-2)=-cost
211       rt(1,2,i-2)=-sint
212       rt(1,3,i-2)=0.0D0
213       rt(2,1,i-2)=sint*cosphi
214       rt(2,2,i-2)=-cost*cosphi
215       rt(2,3,i-2)=sinphi
216       rt(3,1,i-2)=-sint*sinphi
217       rt(3,2,i-2)=cost*sinphi
218       rt(3,3,i-2)=cosphi
219       call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
220       do j=1,3
221         dc_norm(j,i-1)=prod(j,1,i-1)
222         dc(j,i-1)=vbld(i)*prod(j,1,i-1)
223         c(j,i)=c(j,i-1)+dc(j,i-1)
224       enddo
225 cd    print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
226
227 C Now calculate the coordinates of SC(i-1)
228 C
229       call locate_side_chain(i-1)
230       return
231       end
232 c-----------------------------------------------------------------------------
233       subroutine locate_side_chain(i)
234
235 C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
236 C
237       implicit real*8 (a-h,o-z)
238       include 'DIMENSIONS'
239       include 'COMMON.CHAIN'
240       include 'COMMON.LOCAL'
241       include 'COMMON.GEO'
242       include 'COMMON.VAR'
243       include 'COMMON.IOUNITS'
244       include 'COMMON.NAMES'
245       include 'COMMON.INTERACT'
246       dimension xx(3)
247
248 c      dsci=dsc(itype(i))
249 c      dsci_inv=dsc_inv(itype(i))
250       dsci=vbld(i+nres)
251       dsci_inv=vbld_inv(i+nres)
252 #ifdef OSF
253       alphi=alph(i)
254       omegi=omeg(i)
255       if (alphi.ne.alphi) alphi=100.0
256       if (omegi.ne.omegi) omegi=-100.0
257 #else
258       alphi=alph(i)
259       omegi=omeg(i)
260 #endif
261       cosalphi=dcos(alphi)
262       sinalphi=dsin(alphi)
263       cosomegi=dcos(omegi)
264       sinomegi=dsin(omegi) 
265       xp= dsci*cosalphi
266       yp= dsci*sinalphi*cosomegi
267       zp=-dsci*sinalphi*sinomegi
268 * Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
269 * X-axis aligned with the vector DC(*,i)
270       theta2=pi-0.5D0*theta(i+1)
271       cost2=dcos(theta2)
272       sint2=dsin(theta2)
273       xx(1)= xp*cost2+yp*sint2
274       xx(2)=-xp*sint2+yp*cost2
275       xx(3)= zp
276 cd    print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
277 cd   &   xp,yp,zp,(xx(k),k=1,3)
278       do j=1,3
279         xloc(j,i)=xx(j)
280       enddo
281 * Bring the SC vectors to the common coordinate system.
282       xx(1)=xloc(1,i)
283       xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
284       xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
285       do j=1,3
286         xrot(j,i)=xx(j)
287       enddo
288       do j=1,3
289         rj=0.0D0
290         do k=1,3
291           rj=rj+prod(j,k,i-1)*xx(k)
292         enddo
293         dc(j,nres+i)=rj
294         dc_norm(j,nres+i)=rj*dsci_inv
295         c(j,nres+i)=c(j,i)+rj
296       enddo
297       return
298       end
299 c------------------------------------------
300       subroutine returnbox
301       include 'DIMENSIONS'
302 #ifdef MPI
303       include 'mpif.h'
304 #endif
305       include 'COMMON.CONTROL'
306       include 'COMMON.VAR'
307       include 'COMMON.MD'
308 #ifndef LANG0
309       include 'COMMON.LANGEVIN'
310 #else
311       include 'COMMON.LANGEVIN.lang0'
312 #endif
313       include 'COMMON.CHAIN'
314       include 'COMMON.DERIV'
315       include 'COMMON.GEO'
316       include 'COMMON.LOCAL'
317       include 'COMMON.INTERACT'
318       include 'COMMON.IOUNITS'
319       include 'COMMON.NAMES'
320       include 'COMMON.TIME1'
321       include 'COMMON.REMD'
322       include 'COMMON.SETUP'
323       include 'COMMON.MUCA'
324       include 'COMMON.HAIRPIN'
325 C change suggested by Ana - begin
326       integer allareout
327 C change suggested by Ana - end
328         j=1
329         chain_beg=1
330 C        do i=1,nres
331 C       write(*,*) 'initial', i,j,c(j,i)
332 C        enddo
333 C change suggested by Ana - begin
334         allareout=1
335 C change suggested by Ana -end
336         do i=1,nres-1
337          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
338           chain_end=i
339           if (allareout.eq.1) then
340             ireturnval=int(c(j,i)/boxxsize)
341             if (c(j,i).le.0) ireturnval=ireturnval-1
342             do k=chain_beg,chain_end
343               c(j,k)=c(j,k)-ireturnval*boxxsize
344               c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
345             enddo
346 C Suggested by Ana
347             if (chain_beg.eq.1) 
348      &      dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
349 C Suggested by Ana -end
350            endif
351            chain_beg=i+1
352            allareout=1
353          else
354           if (int(c(j,i)/boxxsize).eq.0) allareout=0
355          endif
356         enddo
357          if (allareout.eq.1) then
358             ireturnval=int(c(j,i)/boxxsize)
359             if (c(j,i).le.0) ireturnval=ireturnval-1
360             do k=chain_beg,nres
361               c(j,k)=c(j,k)-ireturnval*boxxsize
362               c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
363             enddo
364           endif
365 C NO JUMP 
366 C        do i=1,nres
367 C        write(*,*) 'befor no jump', i,j,c(j,i)
368 C        enddo
369         nojumpval=0
370         do i=2,nres
371            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
372              difference=abs(c(j,i-1)-c(j,i))
373 C             print *,'diff', difference
374              if (difference.gt.boxxsize/2.0) then
375                 if (c(j,i-1).gt.c(j,i)) then
376                   nojumpval=1
377                  else
378                    nojumpval=-1
379                  endif
380               else
381               nojumpval=0
382               endif
383               endif
384               c(j,i)=c(j,i)+nojumpval*boxxsize
385               c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
386          enddo
387        nojumpval=0
388         do i=2,nres
389            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
390              difference=abs(c(j,i-1)-c(j,i))
391              if (difference.gt.boxxsize/2.0) then
392                 if (c(j,i-1).gt.c(j,i)) then
393                   nojumpval=1
394                  else
395                    nojumpval=-1
396                  endif
397               else
398               nojumpval=0
399               endif
400              endif
401               c(j,i)=c(j,i)+nojumpval*boxxsize
402               c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
403          enddo
404
405 C        do i=1,nres
406 C        write(*,*) 'after no jump', i,j,c(j,i)
407 C        enddo
408
409 C NOW Y dimension
410 C suggesed by Ana begins
411         allareout=1
412 C suggested by Ana ends
413         j=2
414         chain_beg=1
415         do i=1,nres-1
416          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
417           chain_end=i
418           if (allareout.eq.1) then
419             ireturnval=int(c(j,i)/boxysize)
420             if (c(j,i).le.0) ireturnval=ireturnval-1
421             do k=chain_beg,chain_end
422               c(j,k)=c(j,k)-ireturnval*boxysize
423              c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
424             enddo
425 C Suggested by Ana
426             if (chain_beg.eq.1)
427      &      dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
428 C Suggested by Ana -end
429            endif
430            chain_beg=i+1
431            allareout=1
432          else
433           if (int(c(j,i)/boxysize).eq.0) allareout=0
434          endif
435         enddo
436          if (allareout.eq.1) then
437             ireturnval=int(c(j,i)/boxysize)
438             if (c(j,i).le.0) ireturnval=ireturnval-1
439             do k=chain_beg,nres
440               c(j,k)=c(j,k)-ireturnval*boxysize
441               c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
442             enddo
443           endif
444         nojumpval=0
445         do i=2,nres
446            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
447              difference=abs(c(j,i-1)-c(j,i))
448              if (difference.gt.boxysize/2.0) then
449                 if (c(j,i-1).gt.c(j,i)) then
450                   nojumpval=1
451                  else
452                    nojumpval=-1
453                  endif
454              else
455               nojumpval=0
456               endif
457            endif
458               c(j,i)=c(j,i)+nojumpval*boxysize
459               c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
460          enddo
461       nojumpval=0
462         do i=2,nres
463            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
464              difference=abs(c(j,i-1)-c(j,i))
465              if (difference.gt.boxysize/2.0) then
466                 if (c(j,i-1).gt.c(j,i)) then
467                   nojumpval=1
468                  else
469                    nojumpval=-1
470                  endif
471               else
472               nojumpval=0
473               endif
474             endif
475               c(j,i)=c(j,i)+nojumpval*boxysize
476               c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
477          enddo
478 C Now Z dimension
479 C Suggested by Ana -begins
480         allareout=1
481 C Suggested by Ana -ends
482        j=3
483         chain_beg=1
484         do i=1,nres-1
485          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
486           chain_end=i
487           if (allareout.eq.1) then
488             ireturnval=int(c(j,i)/boxysize)
489             if (c(j,i).le.0) ireturnval=ireturnval-1
490             do k=chain_beg,chain_end
491               c(j,k)=c(j,k)-ireturnval*boxzsize
492               c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
493             enddo
494 C Suggested by Ana
495             if (chain_beg.eq.1)
496      &      dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
497 C Suggested by Ana -end
498            endif
499            chain_beg=i+1
500            allareout=1
501          else
502           if (int(c(j,i)/boxzsize).eq.0) allareout=0
503          endif
504         enddo
505          if (allareout.eq.1) then
506             ireturnval=int(c(j,i)/boxzsize)
507             if (c(j,i).le.0) ireturnval=ireturnval-1
508             do k=chain_beg,nres
509               c(j,k)=c(j,k)-ireturnval*boxzsize
510               c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
511             enddo
512           endif
513         nojumpval=0
514         do i=2,nres
515            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
516              difference=abs(c(j,i-1)-c(j,i))
517              if (difference.gt.(boxzsize/2.0)) then
518                 if (c(j,i-1).gt.c(j,i)) then
519                   nojumpval=1
520                  else
521                    nojumpval=-1
522                  endif
523               else
524               nojumpval=0
525               endif
526             endif
527               c(j,i)=c(j,i)+nojumpval*boxzsize
528               c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
529          enddo
530        nojumpval=0
531         do i=2,nres
532            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
533              difference=abs(c(j,i-1)-c(j,i))
534              if (difference.gt.boxzsize/2.0) then
535                 if (c(j,i-1).gt.c(j,i)) then
536                   nojumpval=1
537                  else
538                    nojumpval=-1
539                  endif
540               else
541               nojumpval=0
542               endif
543             endif
544              c(j,i)=c(j,i)+nojumpval*boxzsize
545               c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
546          enddo
547
548         return
549         end
550