3 C Build the virtual polypeptide chain. Side-chain centroids are moveable.
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
26 C first three CAs and SC(2).
28 subroutine chainbuild_extconf
30 C Build the virtual polypeptide chain. Side-chain centroids are moveable.
35 include 'COMMON.CHAIN'
36 include 'COMMON.LOCAL'
39 include 'COMMON.IOUNITS'
40 include 'COMMON.NAMES'
41 include 'COMMON.INTERACT'
43 double precision e1(3),e2(3),e3(3)
44 double precision be,be1,alfai
45 double precision xp,yp,zp,cost2,sint2,cosomegi,sinomegi
46 double precision dist,alpha,beta
47 logical lprn,perbox,fail
50 c write (iout,*) "Calling chainbuild_extconf"
53 * Build the alpha-carbon chain.
56 call locate_next_res(i)
59 C First and last SC must coincide with the corresponding CA.
63 dc_norm(j,nres+1)=0.0D0
65 dc_norm(j,nres+nres)=0.0D0
67 c(j,nres+nres)=c(j,nres)
75 write (iout,'(/a)') 'Recalculated internal coordinates'
78 c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
81 if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i)
82 be1=rad2deg*beta(nres+i,i,maxres2,i+1)
84 if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
85 write (iout,1212) restyp(itype(i)),i,dist(i-1,i),
86 & alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,maxres2),be1
88 1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
95 c-------------------------------------------------------------------------
98 C Define the origin and orientation of the coordinate system and locate
99 C the first three atoms.
103 double precision cost,sint
105 include 'COMMON.CHAIN'
106 include 'COMMON.LOCAL'
160 dc_norm(j,2)=prod(j,1,2)
161 dc(j,2)=vbld(3)*prod(j,1,2)
162 c(j,3)=c(j,2)+dc(j,2)
164 call locate_side_chain(2)
167 c-----------------------------------------------------------------------------
168 subroutine locate_next_res(i)
170 C Locate CA(i) and SC(i-1)
174 double precision theti,phii,cost,sint,cosphi,sinphi
176 include 'COMMON.CHAIN'
177 include 'COMMON.LOCAL'
180 include 'COMMON.IOUNITS'
181 include 'COMMON.NAMES'
182 include 'COMMON.INTERACT'
184 C Define the rotation matrices corresponding to CA(i)
188 if (theti.ne.theti) theti=100.0
190 if (phii.ne.phii) phii=180.0
199 c write (iout,*) "locate_next_res i",i
200 * Define the matrices of the rotation about the virtual-bond valence angles
201 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
202 * program), R(i,j,k), and, the cumulative matrices of rotation RT
224 rt(2,1,i-2)=-sint*cosphi
225 rt(2,2,i-2)=-cost*cosphi
227 rt(3,1,i-2)=-sint*sinphi
228 rt(3,2,i-2)=-cost*sinphi
230 call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
231 c write (iout,*) "prod",i-2
233 c write (iout,*) (prod(j,k,i-2),k=1,3)
235 c write (iout,*) "prod",i-1
237 c write (iout,*) (prod(j,k,i-1),k=1,3)
240 dc_norm(j,i-1)=prod(j,1,i-1)
241 dc(j,i-1)=vbld(i)*prod(j,1,i-1)
242 c(j,i)=c(j,i-1)+dc(j,i-1)
244 c write (iout,*) "dc",i-1,(dc(j,i-1),j=1,3)
245 c write (iout,*) "c",i,(dc(j,i),j=1,3)
246 cd print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
248 C Now calculate the coordinates of SC(i-1)
250 call locate_side_chain(i-1)
253 c-----------------------------------------------------------------------------
254 subroutine sc_coord_rebuild(i)
256 C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
260 include 'COMMON.CHAIN'
261 include 'COMMON.LOCAL'
264 include 'COMMON.IOUNITS'
265 include 'COMMON.NAMES'
266 include 'COMMON.INTERACT'
268 double precision xx(3)
269 double precision dsci,dsci_inv,alphi,omegi,cosalphi,sinalphi,
270 & cosomegi,sinomegi,rp(3),theta2,cost2,sint2,rj
271 double precision scalar
272 double precision ref(3,3),scalp,sscalp,refnorm
273 double precision alpha,beta
275 c dsci_inv=dsc_inv(itype(i))
277 dsci_inv=vbld_inv(i+nres)
281 if (alphi.ne.alphi) alphi=100.0
282 if (omegi.ne.omegi) omegi=-100.0
292 rp(2)= sinalphi*cosomegi
293 rp(3)=-sinalphi*sinomegi
294 c Build the reference system
296 ref(j,1)=-dc_norm(j,i-1)+dc_norm(j,i)
298 refnorm=dsqrt(scalar(ref(1,1),ref(1,1)))
300 ref(j,1)=ref(j,1)/refnorm
302 scalp=scalar(ref(1,1),dc_norm(1,i))
303 sscalp=1.0d0/dsqrt(1.0d0-scalp*scalp)
305 ref(j,2)=(dc_norm(j,i)-scalp*ref(j,1))*sscalp
307 ref(1,3)= ref(2,1)*ref(3,2)-ref(3,1)*ref(2,2)
308 ref(2,3)=-ref(1,1)*ref(3,2)+ref(3,1)*ref(1,2)
309 ref(3,3)= ref(1,1)*ref(2,2)-ref(2,1)*ref(1,2)
311 c write (iout,*) j,scalar(ref(1,j),ref(1,1)),
312 c & scalar(ref(1,j),ref(1,2)),scalar(ref(1,j),ref(1,3))
314 c Bring the coordinates to the global reference system
316 dc_norm(j,nres+i)=0.0d0
318 dc_norm(j,nres+i)=dc_norm(j,nres+i)+ref(j,k)*rp(k)
320 dc(j,nres+i)=dc_norm(j,nres+i)*dsci
321 c(j,nres+i)=c(j,i)+dc(j,nres+i)
323 c write (iout,*) scalar(dc_norm(1,i+nres),dc_norm(1,i+nres)),
324 c & dsqrt(scalar(dc(1,i+nres),dc(1,i+nres)))
325 c Check the internal coordinates
326 c c(:,2*nres+1)=ref(:,1)+c(:,i)
327 c write (iout,*) "alpha",rad2deg*alphi,
328 c & rad2deg*alpha(nres+i,i,2*nres+1)
329 c write (iout,*) "omega",rad2deg*omegi,
330 c & rad2deg*beta(nres+i,i,2*nres+1,i+1)
333 c-----------------------------------------------------------------------------
334 subroutine locate_side_chain(i)
336 C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
340 include 'COMMON.CHAIN'
341 include 'COMMON.LOCAL'
344 include 'COMMON.IOUNITS'
345 include 'COMMON.NAMES'
346 include 'COMMON.INTERACT'
348 double precision xx(3)
349 double precision dsci,dsci_inv,alphi,omegi,cosalphi,sinalphi,
350 & cosomegi,sinomegi,xp,yp,zp,theta2,cost2,sint2,rj
353 c dsci_inv=dsc_inv(itype(i))
355 dsci_inv=vbld_inv(i+nres)
359 if (alphi.ne.alphi) alphi=100.0
360 if (omegi.ne.omegi) omegi=-100.0
370 yp= dsci*sinalphi*cosomegi
371 zp=-dsci*sinalphi*sinomegi
372 * Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
373 * X-axis aligned with the vector DC(*,i)
374 theta2=pi-0.5D0*theta(i+1)
377 xx(1)= xp*cost2+yp*sint2
378 xx(2)=-xp*sint2+yp*cost2
380 cd print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
381 cd & xp,yp,zp,(xx(k),k=1,3)
385 * Bring the SC vectors to the common coordinate system.
387 xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
388 xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
395 rj=rj+prod(j,k,i-1)*xx(k)
398 dc_norm(j,nres+i)=rj*dsci_inv
399 c(j,nres+i)=c(j,i)+rj
403 c------------------------------------------
409 include 'COMMON.CONTROL'
413 include 'COMMON.LANGEVIN'
416 include 'COMMON.LANGEVIN.lang0.5diag'
418 include 'COMMON.LANGEVIN.lang0'
421 include 'COMMON.CHAIN'
422 include 'COMMON.DERIV'
424 include 'COMMON.LOCAL'
425 include 'COMMON.INTERACT'
426 include 'COMMON.IOUNITS'
427 include 'COMMON.NAMES'
428 include 'COMMON.TIME1'
429 include 'COMMON.REMD'
430 include 'COMMON.SETUP'
431 include 'COMMON.MUCA'
432 include 'COMMON.HAIRPIN'
433 C change suggested by Ana - begin
436 C change suggested by Ana - end
440 C write(*,*) 'initial', i,j,c(j,i)
442 C change suggested by Ana - begin
444 C change suggested by Ana -end
446 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
448 if (allareout.eq.1) then
449 ireturnval=int(c(j,i)/boxxsize)
450 if (c(j,i).le.0) ireturnval=ireturnval-1
451 do k=chain_beg,chain_end
452 c(j,k)=c(j,k)-ireturnval*boxxsize
453 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
457 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
458 C Suggested by Ana -end
463 if (int(c(j,i)/boxxsize).eq.0) allareout=0
466 if (allareout.eq.1) then
467 ireturnval=int(c(j,i)/boxxsize)
468 if (c(j,i).le.0) ireturnval=ireturnval-1
470 c(j,k)=c(j,k)-ireturnval*boxxsize
471 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
476 C write(*,*) 'befor no jump', i,j,c(j,i)
480 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
481 difference=abs(c(j,i-1)-c(j,i))
482 C print *,'diff', difference
483 if (difference.gt.boxxsize/2.0) then
484 if (c(j,i-1).gt.c(j,i)) then
493 c(j,i)=c(j,i)+nojumpval*boxxsize
494 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
498 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
499 difference=abs(c(j,i-1)-c(j,i))
500 if (difference.gt.boxxsize/2.0) then
501 if (c(j,i-1).gt.c(j,i)) then
510 c(j,i)=c(j,i)+nojumpval*boxxsize
511 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
515 C write(*,*) 'after no jump', i,j,c(j,i)
519 C suggesed by Ana begins
521 C suggested by Ana ends
525 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
527 if (allareout.eq.1) then
528 ireturnval=int(c(j,i)/boxysize)
529 if (c(j,i).le.0) ireturnval=ireturnval-1
530 do k=chain_beg,chain_end
531 c(j,k)=c(j,k)-ireturnval*boxysize
532 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
536 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
537 C Suggested by Ana -end
542 if (int(c(j,i)/boxysize).eq.0) allareout=0
545 if (allareout.eq.1) then
546 ireturnval=int(c(j,i)/boxysize)
547 if (c(j,i).le.0) ireturnval=ireturnval-1
549 c(j,k)=c(j,k)-ireturnval*boxysize
550 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
555 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
556 difference=abs(c(j,i-1)-c(j,i))
557 if (difference.gt.boxysize/2.0) then
558 if (c(j,i-1).gt.c(j,i)) then
567 c(j,i)=c(j,i)+nojumpval*boxysize
568 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
572 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
573 difference=abs(c(j,i-1)-c(j,i))
574 if (difference.gt.boxysize/2.0) then
575 if (c(j,i-1).gt.c(j,i)) then
584 c(j,i)=c(j,i)+nojumpval*boxysize
585 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
588 C Suggested by Ana -begins
590 C Suggested by Ana -ends
594 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
596 if (allareout.eq.1) then
597 ireturnval=int(c(j,i)/boxysize)
598 if (c(j,i).le.0) ireturnval=ireturnval-1
599 do k=chain_beg,chain_end
600 c(j,k)=c(j,k)-ireturnval*boxzsize
601 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
605 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
606 C Suggested by Ana -end
611 if (int(c(j,i)/boxzsize).eq.0) allareout=0
614 if (allareout.eq.1) then
615 ireturnval=int(c(j,i)/boxzsize)
616 if (c(j,i).le.0) ireturnval=ireturnval-1
618 c(j,k)=c(j,k)-ireturnval*boxzsize
619 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
624 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
625 difference=abs(c(j,i-1)-c(j,i))
626 if (difference.gt.(boxzsize/2.0)) then
627 if (c(j,i-1).gt.c(j,i)) then
636 c(j,i)=c(j,i)+nojumpval*boxzsize
637 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
641 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
642 difference=abs(c(j,i-1)-c(j,i))
643 if (difference.gt.boxzsize/2.0) then
644 if (c(j,i-1).gt.c(j,i)) then
653 c(j,i)=c(j,i)+nojumpval*boxzsize
654 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize