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 locate_side_chain(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,xp,yp,zp,theta2,cost2,sint2,rj
262 c dsci_inv=dsc_inv(itype(i))
264 dsci_inv=vbld_inv(i+nres)
268 if (alphi.ne.alphi) alphi=100.0
269 if (omegi.ne.omegi) omegi=-100.0
279 yp= dsci*sinalphi*cosomegi
280 zp=-dsci*sinalphi*sinomegi
281 * Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
282 * X-axis aligned with the vector DC(*,i)
283 theta2=pi-0.5D0*theta(i+1)
286 xx(1)= xp*cost2+yp*sint2
287 xx(2)=-xp*sint2+yp*cost2
289 cd print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
290 cd & xp,yp,zp,(xx(k),k=1,3)
294 * Bring the SC vectors to the common coordinate system.
296 xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
297 xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
304 rj=rj+prod(j,k,i-1)*xx(k)
307 dc_norm(j,nres+i)=rj*dsci_inv
308 c(j,nres+i)=c(j,i)+rj
312 c------------------------------------------
318 include 'COMMON.CONTROL'
322 include 'COMMON.LANGEVIN'
325 include 'COMMON.LANGEVIN.lang0.5diag'
327 include 'COMMON.LANGEVIN.lang0'
330 include 'COMMON.CHAIN'
331 include 'COMMON.DERIV'
333 include 'COMMON.LOCAL'
334 include 'COMMON.INTERACT'
335 include 'COMMON.IOUNITS'
336 include 'COMMON.NAMES'
337 include 'COMMON.TIME1'
338 include 'COMMON.REMD'
339 include 'COMMON.SETUP'
340 include 'COMMON.MUCA'
341 include 'COMMON.HAIRPIN'
342 C change suggested by Ana - begin
345 C change suggested by Ana - end
349 C write(*,*) 'initial', i,j,c(j,i)
351 C change suggested by Ana - begin
353 C change suggested by Ana -end
355 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
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,chain_end
361 c(j,k)=c(j,k)-ireturnval*boxxsize
362 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
366 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
367 C Suggested by Ana -end
372 if (int(c(j,i)/boxxsize).eq.0) allareout=0
375 if (allareout.eq.1) then
376 ireturnval=int(c(j,i)/boxxsize)
377 if (c(j,i).le.0) ireturnval=ireturnval-1
379 c(j,k)=c(j,k)-ireturnval*boxxsize
380 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
385 C write(*,*) 'befor no jump', i,j,c(j,i)
389 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
390 difference=abs(c(j,i-1)-c(j,i))
391 C print *,'diff', difference
392 if (difference.gt.boxxsize/2.0) then
393 if (c(j,i-1).gt.c(j,i)) then
402 c(j,i)=c(j,i)+nojumpval*boxxsize
403 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
407 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
408 difference=abs(c(j,i-1)-c(j,i))
409 if (difference.gt.boxxsize/2.0) then
410 if (c(j,i-1).gt.c(j,i)) then
419 c(j,i)=c(j,i)+nojumpval*boxxsize
420 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
424 C write(*,*) 'after no jump', i,j,c(j,i)
428 C suggesed by Ana begins
430 C suggested by Ana ends
434 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
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,chain_end
440 c(j,k)=c(j,k)-ireturnval*boxysize
441 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
445 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
446 C Suggested by Ana -end
451 if (int(c(j,i)/boxysize).eq.0) allareout=0
454 if (allareout.eq.1) then
455 ireturnval=int(c(j,i)/boxysize)
456 if (c(j,i).le.0) ireturnval=ireturnval-1
458 c(j,k)=c(j,k)-ireturnval*boxysize
459 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
464 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
465 difference=abs(c(j,i-1)-c(j,i))
466 if (difference.gt.boxysize/2.0) then
467 if (c(j,i-1).gt.c(j,i)) then
476 c(j,i)=c(j,i)+nojumpval*boxysize
477 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
481 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
482 difference=abs(c(j,i-1)-c(j,i))
483 if (difference.gt.boxysize/2.0) then
484 if (c(j,i-1).gt.c(j,i)) then
493 c(j,i)=c(j,i)+nojumpval*boxysize
494 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
497 C Suggested by Ana -begins
499 C Suggested by Ana -ends
503 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
505 if (allareout.eq.1) then
506 ireturnval=int(c(j,i)/boxysize)
507 if (c(j,i).le.0) ireturnval=ireturnval-1
508 do k=chain_beg,chain_end
509 c(j,k)=c(j,k)-ireturnval*boxzsize
510 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
514 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
515 C Suggested by Ana -end
520 if (int(c(j,i)/boxzsize).eq.0) allareout=0
523 if (allareout.eq.1) then
524 ireturnval=int(c(j,i)/boxzsize)
525 if (c(j,i).le.0) ireturnval=ireturnval-1
527 c(j,k)=c(j,k)-ireturnval*boxzsize
528 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
533 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
534 difference=abs(c(j,i-1)-c(j,i))
535 if (difference.gt.(boxzsize/2.0)) then
536 if (c(j,i-1).gt.c(j,i)) then
545 c(j,i)=c(j,i)+nojumpval*boxzsize
546 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
550 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
551 difference=abs(c(j,i-1)-c(j,i))
552 if (difference.gt.boxzsize/2.0) then
553 if (c(j,i-1).gt.c(j,i)) then
562 c(j,i)=c(j,i)+nojumpval*boxzsize
563 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize