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 * Define the matrices of the rotation about the virtual-bond valence angles
200 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
201 * program), R(i,j,k), and, the cumulative matrices of rotation RT
223 rt(2,1,i-2)=sint*cosphi
224 rt(2,2,i-2)=-cost*cosphi
226 rt(3,1,i-2)=-sint*sinphi
227 rt(3,2,i-2)=cost*sinphi
229 call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
231 dc_norm(j,i-1)=prod(j,1,i-1)
232 dc(j,i-1)=vbld(i)*prod(j,1,i-1)
233 c(j,i)=c(j,i-1)+dc(j,i-1)
235 cd print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
237 C Now calculate the coordinates of SC(i-1)
239 call locate_side_chain(i-1)
242 c-----------------------------------------------------------------------------
243 subroutine sc_coord_rebuild(i)
245 C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
249 include 'COMMON.CHAIN'
250 include 'COMMON.LOCAL'
253 include 'COMMON.IOUNITS'
254 include 'COMMON.NAMES'
255 include 'COMMON.INTERACT'
257 double precision xx(3)
258 double precision dsci,dsci_inv,alphi,omegi,cosalphi,sinalphi,
259 & cosomegi,sinomegi,rp(3),theta2,cost2,sint2,rj
260 double precision scalar
261 double precision ref(3,3),scalp,sscalp,refnorm
262 double precision alpha,beta
264 c dsci_inv=dsc_inv(itype(i))
266 dsci_inv=vbld_inv(i+nres)
270 if (alphi.ne.alphi) alphi=100.0
271 if (omegi.ne.omegi) omegi=-100.0
281 rp(2)= sinalphi*cosomegi
282 rp(3)=-sinalphi*sinomegi
283 c Build the reference system
285 ref(j,1)=-dc_norm(j,i-1)+dc_norm(j,i)
287 refnorm=dsqrt(scalar(ref(1,1),ref(1,1)))
289 ref(j,1)=ref(j,1)/refnorm
291 scalp=scalar(ref(1,1),dc_norm(1,i))
292 sscalp=1.0d0/dsqrt(1.0d0-scalp*scalp)
294 ref(j,2)=(dc_norm(j,i)-scalp*ref(j,1))*sscalp
296 ref(1,3)= ref(2,1)*ref(3,2)-ref(3,1)*ref(2,2)
297 ref(2,3)=-ref(1,1)*ref(3,2)+ref(3,1)*ref(1,2)
298 ref(3,3)= ref(1,1)*ref(2,2)-ref(2,1)*ref(1,2)
300 c write (iout,*) j,scalar(ref(1,j),ref(1,1)),
301 c & scalar(ref(1,j),ref(1,2)),scalar(ref(1,j),ref(1,3))
303 c Bring the coordinates to the global reference system
305 dc_norm(j,nres+i)=0.0d0
307 dc_norm(j,nres+i)=dc_norm(j,nres+i)+ref(j,k)*rp(k)
309 dc(j,nres+i)=dc_norm(j,nres+i)*dsci
310 c(j,nres+i)=c(j,i)+dc(j,nres+i)
312 c write (iout,*) scalar(dc_norm(1,i+nres),dc_norm(1,i+nres)),
313 c & dsqrt(scalar(dc(1,i+nres),dc(1,i+nres)))
314 c Check the internal coordinates
315 c c(:,2*nres+1)=ref(:,1)+c(:,i)
316 c write (iout,*) "alpha",rad2deg*alphi,
317 c & rad2deg*alpha(nres+i,i,2*nres+1)
318 c write (iout,*) "omega",rad2deg*omegi,
319 c & rad2deg*beta(nres+i,i,2*nres+1,i+1)
322 c-----------------------------------------------------------------------------
323 subroutine locate_side_chain(i)
325 C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
329 include 'COMMON.CHAIN'
330 include 'COMMON.LOCAL'
333 include 'COMMON.IOUNITS'
334 include 'COMMON.NAMES'
335 include 'COMMON.INTERACT'
337 double precision xx(3)
338 double precision dsci,dsci_inv,alphi,omegi,cosalphi,sinalphi,
339 & cosomegi,sinomegi,xp,yp,zp,theta2,cost2,sint2,rj
342 c dsci_inv=dsc_inv(itype(i))
344 dsci_inv=vbld_inv(i+nres)
348 if (alphi.ne.alphi) alphi=100.0
349 if (omegi.ne.omegi) omegi=-100.0
359 yp= dsci*sinalphi*cosomegi
360 zp=-dsci*sinalphi*sinomegi
361 * Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
362 * X-axis aligned with the vector DC(*,i)
363 theta2=pi-0.5D0*theta(i+1)
366 xx(1)= xp*cost2+yp*sint2
367 xx(2)=-xp*sint2+yp*cost2
369 cd print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
370 cd & xp,yp,zp,(xx(k),k=1,3)
374 * Bring the SC vectors to the common coordinate system.
376 xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
377 xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
384 rj=rj+prod(j,k,i-1)*xx(k)
387 dc_norm(j,nres+i)=rj*dsci_inv
388 c(j,nres+i)=c(j,i)+rj
392 c------------------------------------------
398 include 'COMMON.CONTROL'
402 include 'COMMON.LANGEVIN'
405 include 'COMMON.LANGEVIN.lang0.5diag'
407 include 'COMMON.LANGEVIN.lang0'
410 include 'COMMON.CHAIN'
411 include 'COMMON.DERIV'
413 include 'COMMON.LOCAL'
414 include 'COMMON.INTERACT'
415 include 'COMMON.IOUNITS'
416 include 'COMMON.NAMES'
417 include 'COMMON.TIME1'
418 include 'COMMON.REMD'
419 include 'COMMON.SETUP'
420 include 'COMMON.MUCA'
421 include 'COMMON.HAIRPIN'
422 C change suggested by Ana - begin
425 C change suggested by Ana - end
429 C write(*,*) 'initial', i,j,c(j,i)
431 C change suggested by Ana - begin
433 C change suggested by Ana -end
435 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
437 if (allareout.eq.1) then
438 ireturnval=int(c(j,i)/boxxsize)
439 if (c(j,i).le.0) ireturnval=ireturnval-1
440 do k=chain_beg,chain_end
441 c(j,k)=c(j,k)-ireturnval*boxxsize
442 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
446 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
447 C Suggested by Ana -end
452 if (int(c(j,i)/boxxsize).eq.0) allareout=0
455 if (allareout.eq.1) then
456 ireturnval=int(c(j,i)/boxxsize)
457 if (c(j,i).le.0) ireturnval=ireturnval-1
459 c(j,k)=c(j,k)-ireturnval*boxxsize
460 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
465 C write(*,*) 'befor no jump', i,j,c(j,i)
469 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
470 difference=abs(c(j,i-1)-c(j,i))
471 C print *,'diff', difference
472 if (difference.gt.boxxsize/2.0) then
473 if (c(j,i-1).gt.c(j,i)) then
482 c(j,i)=c(j,i)+nojumpval*boxxsize
483 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
487 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
488 difference=abs(c(j,i-1)-c(j,i))
489 if (difference.gt.boxxsize/2.0) then
490 if (c(j,i-1).gt.c(j,i)) then
499 c(j,i)=c(j,i)+nojumpval*boxxsize
500 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
504 C write(*,*) 'after no jump', i,j,c(j,i)
508 C suggesed by Ana begins
510 C suggested by Ana ends
514 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
516 if (allareout.eq.1) then
517 ireturnval=int(c(j,i)/boxysize)
518 if (c(j,i).le.0) ireturnval=ireturnval-1
519 do k=chain_beg,chain_end
520 c(j,k)=c(j,k)-ireturnval*boxysize
521 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
525 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
526 C Suggested by Ana -end
531 if (int(c(j,i)/boxysize).eq.0) allareout=0
534 if (allareout.eq.1) then
535 ireturnval=int(c(j,i)/boxysize)
536 if (c(j,i).le.0) ireturnval=ireturnval-1
538 c(j,k)=c(j,k)-ireturnval*boxysize
539 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
544 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
545 difference=abs(c(j,i-1)-c(j,i))
546 if (difference.gt.boxysize/2.0) then
547 if (c(j,i-1).gt.c(j,i)) then
556 c(j,i)=c(j,i)+nojumpval*boxysize
557 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
561 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
562 difference=abs(c(j,i-1)-c(j,i))
563 if (difference.gt.boxysize/2.0) then
564 if (c(j,i-1).gt.c(j,i)) then
573 c(j,i)=c(j,i)+nojumpval*boxysize
574 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
577 C Suggested by Ana -begins
579 C Suggested by Ana -ends
583 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
585 if (allareout.eq.1) then
586 ireturnval=int(c(j,i)/boxysize)
587 if (c(j,i).le.0) ireturnval=ireturnval-1
588 do k=chain_beg,chain_end
589 c(j,k)=c(j,k)-ireturnval*boxzsize
590 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
594 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
595 C Suggested by Ana -end
600 if (int(c(j,i)/boxzsize).eq.0) allareout=0
603 if (allareout.eq.1) then
604 ireturnval=int(c(j,i)/boxzsize)
605 if (c(j,i).le.0) ireturnval=ireturnval-1
607 c(j,k)=c(j,k)-ireturnval*boxzsize
608 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
613 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
614 difference=abs(c(j,i-1)-c(j,i))
615 if (difference.gt.(boxzsize/2.0)) then
616 if (c(j,i-1).gt.c(j,i)) then
625 c(j,i)=c(j,i)+nojumpval*boxzsize
626 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
630 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
631 difference=abs(c(j,i-1)-c(j,i))
632 if (difference.gt.boxzsize/2.0) then
633 if (c(j,i-1).gt.c(j,i)) then
642 c(j,i)=c(j,i)+nojumpval*boxzsize
643 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize