2 !-----------------------------------------------------------------------------
11 !-----------------------------------------------------------------------------
15 real(kind=8),dimension(:,:,:),allocatable :: t,r !(3,3,maxres)
16 !-----------------------------------------------------------------------------
19 !-----------------------------------------------------------------------------
21 ! Variables (set in init routine) never modified by local_move
23 integer :: init_called
24 logical :: locmove_output
25 real(kind=8) :: min_theta, max_theta
26 real(kind=8) :: dmin2,dmax2
27 real(kind=8) :: flag,small,small2
28 ! Workspace for local_move
30 integer :: a_n,b_n,res_n
31 real(kind=8),dimension(0:7) :: a_ang
32 real(kind=8),dimension(0:3) :: b_ang
33 real(kind=8),dimension(0:11) :: res_ang
34 logical,dimension(0:2,0:7) :: a_tab
35 logical,dimension(0:2,0:3) :: b_tab
36 logical,dimension(0:2,0:2,0:11) :: res_tab
37 !-----------------------------------------------------------------------------
38 ! integer,dimension(:),allocatable :: itype_pdb !(maxres) initialize in molread
39 !-----------------------------------------------------------------------------
42 !-----------------------------------------------------------------------------
44 !-----------------------------------------------------------------------------
46 !-----------------------------------------------------------------------------
47 real(kind=8) function ARCOS(X)
48 ! implicit real*8 (a-h,o-z)
49 ! include 'COMMON.GEO'
52 IF (DABS(X).LT.1.0D0) GOTO 1
53 ARCOS=PIPOL*(1.0d0-DSIGN(1.0D0,X))
58 !-----------------------------------------------------------------------------
60 !-----------------------------------------------------------------------------
63 ! Build the virtual polypeptide chain. Side-chain centroids are moveable.
66 ! implicit real*8 (a-h,o-z)
67 ! include 'DIMENSIONS'
68 ! include 'COMMON.CHAIN'
69 ! include 'COMMON.LOCAL'
70 ! include 'COMMON.GEO'
71 ! include 'COMMON.VAR'
72 ! include 'COMMON.IOUNITS'
73 ! include 'COMMON.NAMES'
74 ! include 'COMMON.INTERACT'
78 real(kind=8) :: be,be1,alfai
81 ! Set lprn=.true. for debugging
84 ! Define the origin and orientation of the coordinate system and locate the
85 ! first three CA's and SC(2).
89 ! Build the alpha-carbon chain.
92 call locate_next_res(i)
95 ! First and last SC must coincide with the corresponding CA.
99 dc_norm(j,nres+1)=0.0D0
100 dc(j,nres+nres)=0.0D0
101 dc_norm(j,nres+nres)=0.0D0
103 c(j,nres+nres)=c(j,nres)
106 ! Temporary diagnosis
111 write (iout,'(/a)') 'Recalculated internal coordinates'
114 c(j,nres2)=0.5D0*(c(j,i-1)+c(j,i+1)) !maxres2=2*maxres
117 if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i)
118 be1=rad2deg*beta(nres+i,i,nres2,i+1)
120 if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
121 write (iout,1212) restyp(itype(i)),i,dist(i-1,i),&
122 alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,nres2),be1
124 1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
129 end subroutine chainbuild
130 !-----------------------------------------------------------------------------
131 subroutine orig_frame
133 ! Define the origin and orientation of the coordinate system and locate
134 ! the first three atoms.
136 ! implicit real*8 (a-h,o-z)
137 ! include 'DIMENSIONS'
138 ! include 'COMMON.CHAIN'
139 ! include 'COMMON.LOCAL'
140 ! include 'COMMON.GEO'
141 ! include 'COMMON.VAR'
144 real(kind=8) :: cost,sint
146 !el allocate(t(3,3,nres)) !(3,3,maxres)
147 !el allocate(r(3,3,nres)) !(3,3,maxres)
148 !el allocate(rt(3,3,nres)) !(3,3,maxres)
149 !el allocate(dc_norm(3,0:2*nres)) !(3,0:maxres2)
150 !el allocate(prod(3,3,nres)) !(3,3,maxres)
203 dc_norm(j,2)=prod(j,1,2)
204 dc(j,2)=vbld(3)*prod(j,1,2)
205 c(j,3)=c(j,2)+dc(j,2)
207 call locate_side_chain(2)
209 end subroutine orig_frame
210 !-----------------------------------------------------------------------------
211 subroutine locate_next_res(i)
213 ! Locate CA(i) and SC(i-1)
215 ! implicit real*8 (a-h,o-z)
216 ! include 'DIMENSIONS'
217 ! include 'COMMON.CHAIN'
218 ! include 'COMMON.LOCAL'
219 ! include 'COMMON.GEO'
220 ! include 'COMMON.VAR'
221 ! include 'COMMON.IOUNITS'
222 ! include 'COMMON.NAMES'
223 ! include 'COMMON.INTERACT'
225 ! Define the rotation matrices corresponding to CA(i)
229 real(kind=8) :: theti,phii
230 real(kind=8) :: cost,sint,cosphi,sinphi
235 call proc_proc(theti,icrc)
236 if(icrc.eq.1)theti=100.0
239 call proc_proc(phii,icrc)
240 if(icrc.eq.1)phii=180.0
243 if (theti.ne.theti) theti=100.0
245 if (phii.ne.phii) phii=180.0
255 ! Define the matrices of the rotation about the virtual-bond valence angles
256 ! theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
257 ! program), R(i,j,k), and, the cumulative matrices of rotation RT
279 rt(2,1,i-2)=sint*cosphi
280 rt(2,2,i-2)=-cost*cosphi
282 rt(3,1,i-2)=-sint*sinphi
283 rt(3,2,i-2)=cost*sinphi
285 call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
287 dc_norm(j,i-1)=prod(j,1,i-1)
288 dc(j,i-1)=vbld(i)*prod(j,1,i-1)
289 c(j,i)=c(j,i-1)+dc(j,i-1)
291 !d print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
293 ! Now calculate the coordinates of SC(i-1)
295 call locate_side_chain(i-1)
297 end subroutine locate_next_res
298 !-----------------------------------------------------------------------------
299 subroutine locate_side_chain(i)
301 ! Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
303 ! implicit real*8 (a-h,o-z)
304 ! include 'DIMENSIONS'
305 ! include 'COMMON.CHAIN'
306 ! include 'COMMON.LOCAL'
307 ! include 'COMMON.GEO'
308 ! include 'COMMON.VAR'
309 ! include 'COMMON.IOUNITS'
310 ! include 'COMMON.NAMES'
311 ! include 'COMMON.INTERACT'
313 real(kind=8),dimension(3) :: xx
314 real(kind=8) :: alphi,omegi,theta2
315 real(kind=8) :: dsci,dsci_inv,sinalphi,cosalphi,cosomegi,sinomegi
316 real(kind=8) :: xp,yp,zp,cost2,sint2,rj
318 ! dsci_inv=dsc_inv(itype(i))
320 dsci_inv=vbld_inv(i+nres)
327 call proc_proc(alphi,icrc)
328 if(icrc.eq.1)alphi=100.0
330 call proc_proc(omegi,icrc)
331 if(icrc.eq.1)omegi=-100.0
333 if (alphi.ne.alphi) alphi=100.0
334 if (omegi.ne.omegi) omegi=-100.0
345 yp= dsci*sinalphi*cosomegi
346 zp=-dsci*sinalphi*sinomegi
347 ! Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
348 ! X-axis aligned with the vector DC(*,i)
349 theta2=pi-0.5D0*theta(i+1)
352 xx(1)= xp*cost2+yp*sint2
353 xx(2)=-xp*sint2+yp*cost2
355 !d print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
356 !d & xp,yp,zp,(xx(k),k=1,3)
360 ! Bring the SC vectors to the common coordinate system.
362 xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
363 xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
370 rj=rj+prod(j,k,i-1)*xx(k)
373 dc_norm(j,nres+i)=rj*dsci_inv
374 c(j,nres+i)=c(j,i)+rj
377 end subroutine locate_side_chain
378 !-----------------------------------------------------------------------------
380 !-----------------------------------------------------------------------------
381 subroutine int_from_cart1(lprn)
382 ! implicit real*8 (a-h,o-z)
383 ! include 'DIMENSIONS'
388 ! include 'COMMON.IOUNITS'
389 ! include 'COMMON.VAR'
390 ! include 'COMMON.CHAIN'
391 ! include 'COMMON.GEO'
392 ! include 'COMMON.INTERACT'
393 ! include 'COMMON.LOCAL'
394 ! include 'COMMON.NAMES'
395 ! include 'COMMON.SETUP'
396 ! include 'COMMON.TIME1'
400 real(kind=8) :: dnorm1,dnorm2,be
403 if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
411 vbld_inv(nres+1)=0.0d0
412 vbld_inv(2*nres)=0.0d0
415 #if defined(PARINT) && defined(MPI)
416 do i=iint_start,iint_end
423 c(j,nres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1 &
424 +(c(j,i+1)-c(j,i))/dnorm2)
428 if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
429 if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
430 tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
432 if (itype(i-1).ne.10) then
433 tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
434 omicron(1,i)=alpha(i-2,i-1,i-1+nres)
435 omicron(2,i)=alpha(i-1+nres,i-1,i)
437 if (itype(i).ne.10) then
438 tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
441 omeg(i)=beta(nres+i,i,nres2,i+1)
442 alph(i)=alpha(nres+i,i,nres2)
443 theta(i+1)=alpha(i-1,i,i+1)
445 vbld_inv(i)=1.0d0/vbld(i)
446 vbld(nres+i)=dist(nres+i,i)
447 if (itype(i).ne.10) then
448 vbld_inv(nres+i)=1.0d0/vbld(nres+i)
450 vbld_inv(nres+i)=0.0d0
453 #if defined(PARINT) && defined(MPI)
454 if (nfgtasks1.gt.1) then
455 !d write(iout,*) "iint_start",iint_start," iint_count",
456 !d & (iint_count(i),i=0,nfgtasks-1)," iint_displ",
457 !d & (iint_displ(i),i=0,nfgtasks-1)
458 !d write (iout,*) "Gather vbld backbone"
461 call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),&
462 MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),&
463 MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
464 !d write (iout,*) "Gather vbld_inv"
466 call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),&
467 MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),&
468 MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
469 !d write (iout,*) "Gather vbld side chain"
471 call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),&
472 MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),&
473 MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
474 !d write (iout,*) "Gather vbld_inv side chain"
476 call MPI_Allgatherv(vbld_inv(iint_start+nres),&
477 iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),&
478 iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
479 !d write (iout,*) "Gather theta"
481 call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),&
482 MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),&
483 MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
484 !d write (iout,*) "Gather phi"
486 call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),&
487 MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),&
488 MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
490 !d write (iout,*) "Gather alph"
492 call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),&
493 MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),&
494 MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
495 !d write (iout,*) "Gather omeg"
497 call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),&
498 MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),&
499 MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
501 time_gather=time_gather+MPI_Wtime()-time00
506 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
511 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
516 write (iout,1212) restyp(itype(i)),i,vbld(i),&
517 rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),&
518 rad2deg*alph(i),rad2deg*omeg(i)
521 1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
523 time_intfcart=time_intfcart+MPI_Wtime()-time01
526 end subroutine int_from_cart1
528 !-----------------------------------------------------------------------------
530 !-----------------------------------------------------------------------------
531 subroutine check_sc_distr
532 ! implicit real*8 (a-h,o-z)
533 ! include 'DIMENSIONS'
534 ! include 'COMMON.TIME1'
535 ! include 'COMMON.INTERACT'
536 ! include 'COMMON.NAMES'
537 ! include 'COMMON.GEO'
538 ! include 'COMMON.HEADER'
539 ! include 'COMMON.CONTROL'
541 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
542 real(kind=8) :: hrtime,mintime,sectime
543 integer,parameter :: MaxSample=10000000
544 real(kind=8),parameter :: delt=1.0D0/MaxSample
545 real(kind=8),dimension(0:72,0:90) :: prob
547 integer :: it,i,j,isample,indal,indom
548 real(kind=8) :: al,om,dV
549 dV=2.0D0*5.0D0*deg2rad*deg2rad
552 if (it.eq.10) goto 10
553 open (20,file=restyp(it)//'_distr.sdc',status='unknown')
554 call gen_side(it,90.0D0 * deg2rad,al,om,fail)
557 open (20,file=restyp(it)//'_distr1.sdc',status='unknown')
563 do isample=1,MaxSample
564 call gen_side(it,90.0D0 * deg2rad,al,om,fail)
566 indom=(rad2deg*om+180.0D0)/5
567 prob(indom,indal)=prob(indom,indal)+delt
571 write (20,'(2f10.3,1pd15.5)') 2*i+0.0D0,5*j-180.0D0,&
577 end subroutine check_sc_distr
579 !-----------------------------------------------------------------------------
581 !-----------------------------------------------------------------------------
582 subroutine geom_to_var(n,x)
584 ! Transfer the geometry parameters to the variable array.
585 ! The positions of variables are as follows:
586 ! 1. Virtual-bond torsional angles: 1 thru nres-3
587 ! 2. Virtual-bond valence angles: nres-2 thru 2*nres-5
588 ! 3. The polar angles alpha of local SC orientation: 2*nres-4 thru
590 ! 4. The torsional angles omega of SC orientation: 2*nres-4+nside+1
591 ! thru 2*nre-4+2*nside
593 ! implicit real*8 (a-h,o-z)
594 ! include 'DIMENSIONS'
595 ! include 'COMMON.VAR'
596 ! include 'COMMON.GEO'
597 ! include 'COMMON.CHAIN'
599 real(kind=8),dimension(n) :: x
600 !d print *,'nres',nres,' nphi',nphi,' ntheta',ntheta,' nvar',nvar
603 !d print *,i,i-3,phi(i)
605 if (n.eq.nphi) return
608 !d print *,i,i-2+nphi,theta(i)
610 if (n.eq.nphi+ntheta) return
612 if (ialph(i,1).gt.0) then
613 x(ialph(i,1))=alph(i)
614 x(ialph(i,1)+nside)=omeg(i)
615 !d print *,i,ialph(i,1),ialph(i,1)+nside,alph(i),omeg(i)
619 end subroutine geom_to_var
620 !-----------------------------------------------------------------------------
621 subroutine var_to_geom(n,x)
623 ! Update geometry parameters according to the variable array.
625 ! implicit real*8 (a-h,o-z)
626 ! include 'DIMENSIONS'
627 ! include 'COMMON.VAR'
628 ! include 'COMMON.CHAIN'
629 ! include 'COMMON.GEO'
630 ! include 'COMMON.IOUNITS'
632 real(kind=8),dimension(n) :: x
633 logical :: change !,reduce
640 if (n.gt.nphi+ntheta) then
643 alph(ii)=x(nphi+ntheta+i)
644 omeg(ii)=pinorm(x(nphi+ntheta+nside+i))
650 if (n.eq.nphi) return
653 if (theta(i).eq.pi) theta(i)=0.99d0*pi
657 end subroutine var_to_geom
658 !-----------------------------------------------------------------------------
659 logical function convert_side(alphi,omegi)
661 real(kind=8) :: alphi,omegi
662 !el real(kind=8) :: pinorm
663 ! include 'COMMON.GEO'
665 ! Apply periodicity restrictions.
666 if (alphi.gt.pi) then
668 omegi=pinorm(omegi+pi)
672 end function convert_side
673 !-----------------------------------------------------------------------------
674 logical function reduce(x)
676 ! Apply periodic restrictions to variables.
678 ! implicit real*8 (a-h,o-z)
679 ! include 'DIMENSIONS'
680 ! include 'COMMON.VAR'
681 ! include 'COMMON.CHAIN'
682 ! include 'COMMON.GEO'
683 logical :: zm,zmiana !,convert_side
684 real(kind=8),dimension(nvar) :: x
688 x(i-3)=pinorm(x(i-3))
690 if (nvar.gt.nphi+ntheta) then
694 x(ii)=thetnorm(x(ii))
695 x(iii)=pinorm(x(iii))
696 ! Apply periodic restrictions.
697 zm=convert_side(x(ii),x(iii))
701 if (nvar.eq.nphi) return
705 x(ii)=dmod(x(ii),dwapi)
706 ! Apply periodic restrictions.
707 if (x(ii).gt.pi) then
710 if (iii.gt.0) x(iii)=pinorm(x(iii)+pi)
711 if (i.lt.nres) x(iii+1)=pinorm(x(iii+1)+pi)
714 x(ii)=dmod(pi-x(ii),dwapi)
715 x(ii+nside)=pinorm(-x(ii+nside))
716 zm=convert_side(x(ii),x(ii+nside))
718 else if (x(ii).lt.-pi) then
723 x(ii)=dmod(pi-x(ii),dwapi)
724 x(ii+nside)=pinorm(-pi-x(ii+nside))
725 zm=convert_side(x(ii),x(ii+nside))
727 else if (x(ii).lt.0.0d0) then
730 if (iii.gt.0) x(iii)=pinorm(x(iii)+pi)
731 if (i.lt.nres) x(iii+1)=pinorm(x(iii+1)+pi)
734 x(ii+nside)=pinorm(-x(ii+nside))
735 zm=convert_side(x(ii),x(ii+nside))
742 !-----------------------------------------------------------------------------
743 real(kind=8) function thetnorm(x)
744 ! This function puts x within [0,2Pi].
747 ! include 'COMMON.GEO'
749 if (xx.lt.0.0d0) xx=xx+dwapi
750 if (xx.gt.0.9999d0*pi) xx=0.9999d0*pi
753 end function thetnorm
755 !-----------------------------------------------------------------------------
756 subroutine var_to_geom_restr(n,xx)
758 ! Update geometry parameters according to the variable array.
760 ! implicit real*8 (a-h,o-z)
761 ! include 'DIMENSIONS'
762 ! include 'COMMON.VAR'
763 ! include 'COMMON.CHAIN'
764 ! include 'COMMON.GEO'
765 ! include 'COMMON.IOUNITS'
767 real(kind=8),dimension(6*nres) :: x,xx !(maxvar) (maxvar=6*maxres)
768 logical :: change !,reduce
774 alph(ii)=x(nphi+ntheta+i)
775 omeg(ii)=pinorm(x(nphi+ntheta+nside+i))
782 if (theta(i).eq.pi) theta(i)=0.99d0*pi
786 end subroutine var_to_geom_restr
787 !-----------------------------------------------------------------------------
789 !-----------------------------------------------------------------------------
790 subroutine gen_rand_conf(nstart,*)
791 ! Generate random conformation or chain cut and regrowth.
793 ! implicit real*8 (a-h,o-z)
794 ! include 'DIMENSIONS'
795 ! include 'COMMON.CHAIN'
796 ! include 'COMMON.LOCAL'
797 ! include 'COMMON.VAR'
798 ! include 'COMMON.INTERACT'
799 ! include 'COMMON.IOUNITS'
800 ! include 'COMMON.MCM'
801 ! include 'COMMON.GEO'
802 ! include 'COMMON.CONTROL'
803 logical :: back,fail !overlap,
805 integer :: i,nstart,maxsi,nsi,maxnit,nit,niter
806 integer :: it1,it2,it,j
807 !d print *,' CG Processor',me,' maxgen=',maxgen
809 !d write (iout,*) 'Gen_Rand_conf: nstart=',nstart
810 if (nstart.lt.5) then
812 phi(4)=gen_phi(4,iabs(itype(2)),iabs(itype(3)))
813 ! write(iout,*)'phi(4)=',rad2deg*phi(4)
814 if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2)),pi,phi(4))
815 ! write(iout,*)'theta(3)=',rad2deg*theta(3)
819 do while (fail.and.nsi.le.maxsi)
820 call gen_side(it1,theta(3),alph(2),omeg(2),fail)
823 if (nsi.gt.maxsi) return 1
838 do while (i.le.nres .and. niter.lt.maxgen)
839 if (i.lt.nstart) then
841 write (iout,'(/80(1h*)/2a/80(1h*))') &
842 'Generation procedure went down to ',&
843 'chain beginning. Cannot continue...'
844 write (*,'(/80(1h*)/2a/80(1h*))') &
845 'Generation procedure went down to ',&
846 'chain beginning. Cannot continue...'
853 ! print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
854 ! & ' nit=',nit,' niter=',niter,' maxgen=',maxgen
855 phi(i+1)=gen_phi(i+1,it1,it)
857 phi(i)=gen_phi(i+1,it2,it1)
858 ! print *,'phi(',i,')=',phi(i)
859 theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
863 do while (fail.and.nsi.le.maxsi)
864 call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail)
867 if (nsi.gt.maxsi) return 1
869 call locate_next_res(i-1)
871 theta(i)=gen_theta(it1,phi(i),phi(i+1))
875 do while (fail.and.nsi.le.maxsi)
876 call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
879 if (nsi.gt.maxsi) return 1
881 call locate_next_res(i)
882 if (overlap(i-1)) then
883 if (nit.lt.maxnit) then
893 'Cannot generate non-overlaping conformation. Increase MAXNIT.'
895 'Cannot generate non-overlaping conformation. Increase MAXNIT.'
906 if (niter.ge.maxgen) then
907 write (iout,'(a,2i5)') &
908 'Too many trials in conformation generation',niter,maxgen
909 write (*,'(a,2i5)') &
910 'Too many trials in conformation generation',niter,maxgen
915 c(j,nres+nres)=c(j,nres)
918 end subroutine gen_rand_conf
919 !-----------------------------------------------------------------------------
920 logical function overlap(i)
921 ! implicit real*8 (a-h,o-z)
922 ! include 'DIMENSIONS'
923 ! include 'COMMON.CHAIN'
924 ! include 'COMMON.INTERACT'
925 ! include 'COMMON.FFIELD'
926 integer :: i,j,iti,itj,iteli,itelj,k
927 real(kind=8) :: redfac,rcomp
933 if (iti.gt.ntyp) return
934 ! Check for SC-SC overlaps.
935 !d print *,'nnt=',nnt,' nct=',nct
938 if (j.lt.i-1 .or. ipot.ne.4) then
939 rcomp=sigmaii(iti,itj)
944 if (dist(nres+i,nres+j).lt.redfac*rcomp) then
946 ! print *,'overlap, SC-SC: i=',i,' j=',j,
947 ! & ' dist=',dist(nres+i,nres+j),' rcomp=',
952 ! Check for overlaps between the added peptide group and the preceding
956 c(j,nres2+1)=0.5D0*(c(j,i)+c(j,i+1))
960 !d print *,'overlap, p-Sc: i=',i,' j=',j,
961 !d & ' dist=',dist(nres+j,maxres2+1)
962 if (dist(nres+j,nres2+1).lt.4.0D0*redfac) then
967 ! Check for overlaps between the added side chain and the preceding peptide
971 c(k,nres2+1)=0.5D0*(c(k,j)+c(k,j+1))
973 !d print *,'overlap, SC-p: i=',i,' j=',j,
974 !d & ' dist=',dist(nres+i,maxres2+1)
975 if (dist(nres+i,nres2+1).lt.4.0D0*redfac) then
980 ! Check for p-p overlaps
982 c(j,nres2+2)=0.5D0*(c(j,i)+c(j,i+1))
987 c(k,nres2+2)=0.5D0*(c(k,j)+c(k,j+1))
989 !d print *,'overlap, p-p: i=',i,' j=',j,
990 !d & ' dist=',dist(maxres2+1,maxres2+2)
991 if(iteli.ne.0.and.itelj.ne.0)then
992 if (dist(nres2+1,nres2+2).lt.rpp(iteli,itelj)*redfac) then
1000 !-----------------------------------------------------------------------------
1001 real(kind=8) function gen_phi(i,it1,it2)
1002 use random, only:ran_number
1003 ! implicit real*8 (a-h,o-z)
1004 ! include 'DIMENSIONS'
1005 ! include 'COMMON.GEO'
1006 ! include 'COMMON.BOUNDS'
1007 integer :: i,it1,it2
1008 ! gen_phi=ran_number(-pi,pi)
1009 ! 8/13/98 Generate phi using pre-defined boundaries
1010 gen_phi=ran_number(phibound(1,i),phibound(2,i))
1012 end function gen_phi
1013 !-----------------------------------------------------------------------------
1014 real(kind=8) function gen_theta(it,gama,gama1)
1015 use random,only:binorm
1016 ! implicit real*8 (a-h,o-z)
1017 ! include 'DIMENSIONS'
1018 ! include 'COMMON.LOCAL'
1019 ! include 'COMMON.GEO'
1020 real(kind=8),dimension(2) :: y,z
1021 real(kind=8) :: theta_max,theta_min,sig,ak
1024 real(kind=8) :: gama,gama1,thet_pred_mean,theta_temp
1025 ! print *,'gen_theta: it=',it
1028 if (dabs(gama).gt.dwapi) then
1035 if (dabs(gama1).gt.dwapi) then
1042 thet_pred_mean=a0thet(it)
1044 thet_pred_mean=thet_pred_mean+athet(k,it,1,1)*y(k) &
1045 +bthet(k,it,1,1)*z(k)
1049 sig=sig*thet_pred_mean+polthet(j,it)
1051 sig=0.5D0/(sig*sig+sigc0(it))
1052 ak=dexp(gthet(1,it)- &
1053 0.5D0*((gthet(2,it)-thet_pred_mean)/gthet(3,it))**2)
1054 ! print '(i5,5(1pe14.4))',it,(gthet(j,it),j=1,3)
1055 ! print '(5(1pe14.4))',thet_pred_mean,theta0(it),sig,sig0(it),ak
1056 theta_temp=binorm(thet_pred_mean,theta0(it),sig,sig0(it),ak)
1057 if (theta_temp.lt.theta_min) theta_temp=theta_min
1058 if (theta_temp.gt.theta_max) theta_temp=theta_max
1059 gen_theta=theta_temp
1060 ! print '(a)','Exiting GENTHETA.'
1062 end function gen_theta
1063 !-----------------------------------------------------------------------------
1064 subroutine gen_side(it,the,al,om,fail)
1065 use random, only:ran_number,mult_norm1
1066 ! implicit real*8 (a-h,o-z)
1067 ! include 'DIMENSIONS'
1068 ! include 'COMMON.GEO'
1069 ! include 'COMMON.LOCAL'
1070 ! include 'COMMON.SETUP'
1071 ! include 'COMMON.IOUNITS'
1072 real(kind=8) :: MaxBoxLen=10.0D0
1073 real(kind=8),dimension(3,3) :: Ap_inv,a,vec
1074 real(kind=8),dimension(:,:),allocatable :: z !(3,maxlob)
1075 real(kind=8),dimension(:),allocatable :: W1,detAp !(maxlob)
1076 real(kind=8),dimension(:),allocatable :: sumW !(0:maxlob)
1077 real(kind=8),dimension(2) :: y,cm,eig
1078 real(kind=8),dimension(2,2) :: box
1079 real(kind=8),dimension(100) :: work
1080 real(kind=8) :: eig_limit=1.0D-8
1081 real(kind=8) :: Big=10.0D0
1082 logical :: lprint,fail,lcheck
1084 integer :: it,i,j,k,l,nlobit,ial,iom,iii,ilob
1085 real(kind=8) :: the,al,om,detApi,wart,y2,wykl,radmax
1086 real(kind=8) :: tant,zz1,W1i,radius,zk,fac,dV,sum,sum1
1087 real(kind=8) :: which_lobe
1091 if (the.eq.0.0D0 .or. the.eq.pi) then
1093 write (*,'(a,i4,a,i3,a,1pe14.5)') &
1094 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the
1096 !d write (iout,'(a,i3,a,1pe14.5)')
1097 !d & 'Error in GenSide: it=',it,' theta=',the
1102 tant=dtan(the-pipol)
1104 allocate(z(3,nlobit))
1105 allocate(W1(nlobit))
1106 allocate(detAp(nlobit))
1107 allocate(sumW(0:nlobit))
1110 print '(a,i4,a)','CG Processor:',me,' Enter Gen_Side.'
1111 write (iout,'(a,i4,a)') 'Processor:',me,' Enter Gen_Side.'
1113 print *,'it=',it,' nlobit=',nlobit,' the=',the,' tant=',tant
1114 write (iout,*) 'it=',it,' nlobit=',nlobit,' the=',the,&
1118 zz1=tant-censc(1,i,it)
1121 a(k,l)=gaussc(k,l,i,it)
1124 detApi=a(2,2)*a(3,3)-a(2,3)**2
1125 Ap_inv(2,2)=a(3,3)/detApi
1126 Ap_inv(2,3)=-a(2,3)/detApi
1127 Ap_inv(3,2)=Ap_inv(2,3)
1128 Ap_inv(3,3)=a(2,2)/detApi
1130 write (*,'(/a,i2/)') 'Cluster #',i
1131 write (*,'(3(1pe14.5),5x,1pe14.5)') &
1132 ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
1133 write (iout,'(/a,i2/)') 'Cluster #',i
1134 write (iout,'(3(1pe14.5),5x,1pe14.5)') &
1135 ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
1140 W1i=W1i+a(k,1)*a(l,1)*Ap_inv(k,l)
1144 W1(i)=dexp(bsc(i,it)-0.5D0*W1i*zz1*zz1)
1145 ! if (lprint) write(*,'(a,3(1pe15.5)/)')
1146 ! & 'detAp, W1, anormi',detApi,W1i,anormi
1150 zk=zk+zz1*Ap_inv(k,l)*a(l,1)
1154 detAp(i)=dsqrt(detApi)
1158 print *,'W1:',(w1(i),i=1,nlobit)
1159 print *,'detAp:',(detAp(i),i=1,nlobit)
1162 print '(i2,3f10.5)',i,(rad2deg*z(j,i),j=2,3)
1164 write (iout,*) 'W1:',(w1(i),i=1,nlobit)
1165 write (iout,*) 'detAp:',(detAp(i),i=1,nlobit)
1168 write (iout,'(i2,3f10.5)') i,(rad2deg*z(j,i),j=2,3)
1172 ! Writing the distribution just to check the procedure
1174 dV=deg2rad**2*10.0D0
1178 fac=fac+W1(i)/detAp(i)
1180 fac=1.0D0/(2.0D0*fac*pi)
1181 !d print *,it,'fac=',fac
1190 a(j-1,k-1)=gaussc(j,k,i,it)
1202 wykl=wykl+a(j,k)*(y(j)-z(j+1,i))*(y(k)-z(k+1,i))
1205 wart=wart+W1(i)*dexp(-0.5D0*wykl)
1212 ! print *,'y',y(1),y(2),' fac=',fac
1214 write (20,'(2f10.3,1pd15.5)') y(1)*rad2deg,y(2)*rad2deg,wart
1219 ! print *,'it=',it,' sum=',sum*dV,' sum1=',sum1*dV
1223 ! Calculate the CM of the system
1226 W1(i)=W1(i)/detAp(i)
1230 sumW(i)=sumW(i-1)+W1(i)
1235 cm(1)=cm(1)+z(2,j)*W1(j)
1236 cm(2)=cm(2)+W1(j)*(z(3,1)+pinorm(z(3,j)-z(3,1)))
1238 cm(1)=cm(1)/sumW(nlobit)
1239 cm(2)=cm(2)/sumW(nlobit)
1240 if (cm(1).gt.Big .or. cm(1).lt.-Big .or. &
1241 cm(2).gt.Big .or. cm(2).lt.-Big) then
1242 !d write (iout,'(a)')
1243 !d & 'Unexpected error in GenSide - CM coordinates too large.'
1244 !d write (iout,'(i5,2(1pe14.5))') it,cm(1),cm(2)
1246 !d & 'Unexpected error in GenSide - CM coordinates too large.'
1247 !d write (*,'(i5,2(1pe14.5))') it,cm(1),cm(2)
1251 !d print *,'CM:',cm(1),cm(2)
1253 ! Find the largest search distance from CM
1259 a(j-1,k-1)=gaussc(j,k,i,it)
1263 call f02faf('N','U',2,a,3,eig,work,100,ifail)
1265 call djacob(2,3,10000,1.0d-10,a,vec,eig)
1269 print *,'*************** CG Processor',me
1270 print *,'CM:',cm(1),cm(2)
1271 write (iout,*) '*************** CG Processor',me
1272 write (iout,*) 'CM:',cm(1),cm(2)
1273 print '(A,8f10.5)','Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
1274 write (iout,'(A,8f10.5)') &
1275 'Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
1278 if (eig(1).lt.eig_limit) then
1280 'From Mult_Norm: Eigenvalues of A are too small.'
1282 'From Mult_Norm: Eigenvalues of A are too small.'
1289 radius=radius+pinorm(z(j+1,i)-cm(j))**2
1291 radius=dsqrt(radius)+3.0D0/dsqrt(eig(1))
1292 if (radius.gt.radmax) radmax=radius
1294 if (radmax.gt.pi) radmax=pi
1296 ! Determine the boundaries of the search rectangle.
1299 print '(a,4(1pe14.4))','W1: ',(W1(i),i=1,nlob(it) )
1300 print '(a,4(1pe14.4))','radmax: ',radmax
1302 box(1,1)=dmax1(cm(1)-radmax,0.0D0)
1303 box(2,1)=dmin1(cm(1)+radmax,pi)
1304 box(1,2)=cm(2)-radmax
1305 box(2,2)=cm(2)+radmax
1308 print *,'CG Processor',me,' Array BOX:'
1310 print *,'Array BOX:'
1312 print '(4(1pe14.4))',((box(k,j),k=1,2),j=1,2)
1313 print '(a,4(1pe14.4))','sumW: ',(sumW(i),i=0,nlob(it) )
1315 write (iout,*)'CG Processor',me,' Array BOX:'
1317 write (iout,*)'Array BOX:'
1319 write(iout,'(4(1pe14.4))') ((box(k,j),k=1,2),j=1,2)
1320 write(iout,'(a,4(1pe14.4))')'sumW: ',(sumW(i),i=0,nlob(it) )
1322 if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
1324 write (iout,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
1325 write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
1327 ! write (iout,'(a)') 'Bad sampling box.'
1332 which_lobe=ran_number(0.0D0,sumW(nlobit))
1333 ! print '(a,1pe14.4)','which_lobe=',which_lobe
1335 if (sumW(i-1).le.which_lobe .and. sumW(i).ge.which_lobe) goto 1
1338 ! print *,'ilob=',ilob,' nlob=',nlob(it)
1342 a(i-1,j-1)=gaussc(i,j,ilob,it)
1345 !d print '(a,i4,a)','CG Processor',me,' Calling MultNorm1.'
1346 call mult_norm1(3,2,a,cm,box,y,fail)
1350 !d print *,'al=',al,' om=',om
1353 end subroutine gen_side
1354 !-----------------------------------------------------------------------------
1355 subroutine overlap_sc(scfail)
1357 ! Internal and cartesian coordinates must be consistent as input,
1358 ! and will be up-to-date on return.
1359 ! At the end of this procedure, scfail is true if there are
1360 ! overlapping residues left, or false otherwise (success)
1362 ! implicit real*8 (a-h,o-z)
1363 ! include 'DIMENSIONS'
1364 ! include 'COMMON.CHAIN'
1365 ! include 'COMMON.INTERACT'
1366 ! include 'COMMON.FFIELD'
1367 ! include 'COMMON.VAR'
1368 ! include 'COMMON.SBRIDGE'
1369 ! include 'COMMON.IOUNITS'
1370 logical :: had_overlaps,fail,scfail
1371 integer,dimension(nres) :: ioverlap !(maxres)
1372 integer :: ioverlap_last,k,maxsi,i,iti,nsi
1375 had_overlaps=.false.
1376 call overlap_sc_list(ioverlap,ioverlap_last)
1377 if (ioverlap_last.gt.0) then
1378 write (iout,*) '#OVERLAPing residues ',ioverlap_last
1379 write (iout,'(20i4)') (ioverlap(k),k=1,ioverlap_last)
1385 if (ioverlap_last.eq.0) exit
1387 do ires=1,ioverlap_last
1393 do while (fail.and.nsi.le.maxsi)
1394 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
1402 call overlap_sc_list(ioverlap,ioverlap_last)
1403 ! write (iout,*) 'Overlaping residues ',ioverlap_last,
1404 ! & (ioverlap(j),j=1,ioverlap_last)
1407 if (k.le.1000.and.ioverlap_last.eq.0) then
1409 if (had_overlaps) then
1410 write (iout,*) '#OVERLAPing all corrected after ',k,&
1411 ' random generation'
1415 write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
1416 write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last)
1422 write (iout,'(a30,i5,a12,i4)') &
1423 '#OVERLAP FAIL in gen_side after',maxsi,&
1427 end subroutine overlap_sc
1428 !-----------------------------------------------------------------------------
1429 subroutine overlap_sc_list(ioverlap,ioverlap_last)
1431 ! implicit real*8 (a-h,o-z)
1432 ! include 'DIMENSIONS'
1433 ! include 'COMMON.GEO'
1434 ! include 'COMMON.LOCAL'
1435 ! include 'COMMON.IOUNITS'
1436 ! include 'COMMON.CHAIN'
1437 ! include 'COMMON.INTERACT'
1438 ! include 'COMMON.FFIELD'
1439 ! include 'COMMON.VAR'
1440 ! include 'COMMON.CALC'
1442 integer,dimension(nres) :: ioverlap !(maxres)
1443 integer :: ioverlap_last
1446 real(kind=8) :: redfac,sig !rrij,sigsq,
1447 integer :: itypi,itypj,itypi1
1448 real(kind=8) :: xi,yi,zi,sig0ij,rcomp,rrij,rij_shift
1452 ! Check for SC-SC overlaps and mark residues
1453 ! print *,'>>overlap_sc nnt=',nnt,' nct=',nct
1455 do i=iatsc_s,iatsc_e
1456 itypi=iabs(itype(i))
1457 itypi1=iabs(itype(i+1))
1461 dxi=dc_norm(1,nres+i)
1462 dyi=dc_norm(2,nres+i)
1463 dzi=dc_norm(3,nres+i)
1464 dsci_inv=dsc_inv(itypi)
1466 do iint=1,nint_gr(i)
1467 do j=istart(i,iint),iend(i,iint)
1469 itypj=iabs(itype(j))
1470 dscj_inv=dsc_inv(itypj)
1471 sig0ij=sigma(itypi,itypj)
1472 chi1=chi(itypi,itypj)
1473 chi2=chi(itypj,itypi)
1480 alf12=0.5D0*(alf1+alf2)
1482 rcomp=sigmaii(itypi,itypj)
1484 rcomp=sigma(itypi,itypj)
1486 ! print '(2(a3,2i3),a3,2f10.5)',
1487 ! & ' i=',i,iti,' j=',j,itj,' d=',dist(nres+i,nres+j)
1492 dxj=dc_norm(1,nres+j)
1493 dyj=dc_norm(2,nres+j)
1494 dzj=dc_norm(3,nres+j)
1495 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1499 sig=sig0ij*dsqrt(sigsq)
1500 rij_shift=1.0D0/rij-sig+sig0ij
1502 !t if ( 1.0/rij .lt. redfac*rcomp .or.
1503 !t & rij_shift.le.0.0D0 ) then
1504 if ( rij_shift.le.0.0D0 ) then
1505 !d write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)')
1506 !d & 'overlap SC-SC: i=',i,' j=',j,
1507 !d & ' dist=',dist(nres+i,nres+j),' rcomp=',
1508 !d & rcomp,1.0/rij,rij_shift
1509 ioverlap_last=ioverlap_last+1
1510 ioverlap(ioverlap_last)=i
1511 do k=1,ioverlap_last-1
1512 if (ioverlap(k).eq.i) ioverlap_last=ioverlap_last-1
1514 ioverlap_last=ioverlap_last+1
1515 ioverlap(ioverlap_last)=j
1516 do k=1,ioverlap_last-1
1517 if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1
1524 end subroutine overlap_sc_list
1526 !-----------------------------------------------------------------------------
1527 ! energy_p_new_barrier.F
1528 !-----------------------------------------------------------------------------
1529 subroutine sc_angular
1530 ! Calculate eps1,eps2,eps3,sigma, and parts of their derivatives in om1,om2,
1531 ! om12. Called by ebp, egb, and egbv.
1534 ! include 'COMMON.CALC'
1535 ! include 'COMMON.IOUNITS'
1539 om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
1540 om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
1541 om12=dxi*dxj+dyi*dyj+dzi*dzj
1543 ! Calculate eps1(om12) and its derivative in om12
1544 faceps1=1.0D0-om12*chiom12
1545 faceps1_inv=1.0D0/faceps1
1546 eps1=dsqrt(faceps1_inv)
1547 ! Following variable is eps1*deps1/dom12
1548 eps1_om12=faceps1_inv*chiom12
1553 ! write (iout,*) "om12",om12," eps1",eps1
1554 ! Calculate sigma(om1,om2,om12) and the derivatives of sigma**2 in om1,om2,
1559 facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
1560 sigsq=1.0D0-facsig*faceps1_inv
1561 sigsq_om1=(chiom1-chiom12*om2)*faceps1_inv
1562 sigsq_om2=(chiom2-chiom12*om1)*faceps1_inv
1563 sigsq_om12=-chi12*(om1om2*faceps1-om12*facsig)*faceps1_inv**2
1569 ! write (iout,*) "chiom1",chiom1," chiom2",chiom2," chiom12",chiom12
1570 ! write (iout,*) "faceps1",faceps1," faceps1_inv",faceps1_inv,
1572 ! Calculate eps2 and its derivatives in om1, om2, and om12.
1575 chipom12=chip12*om12
1576 facp=1.0D0-om12*chipom12
1578 facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12
1579 ! write (iout,*) "chipom1",chipom1," chipom2",chipom2,
1580 ! & " chipom12",chipom12," facp",facp," facp_inv",facp_inv
1581 ! Following variable is the square root of eps2
1582 eps2rt=1.0D0-facp1*facp_inv
1583 ! Following three variables are the derivatives of the square root of eps
1584 ! in om1, om2, and om12.
1585 eps2rt_om1=-4.0D0*(chipom1-chipom12*om2)*facp_inv
1586 eps2rt_om2=-4.0D0*(chipom2-chipom12*om1)*facp_inv
1587 eps2rt_om12=4.0D0*chip12*(om1om2*facp-om12*facp1)*facp_inv**2
1588 ! Evaluate the "asymmetric" factor in the VDW constant, eps3
1589 eps3rt=1.0D0-alf1*om1+alf2*om2-alf12*om12
1590 ! write (iout,*) "eps2rt",eps2rt," eps3rt",eps3rt
1591 ! write (iout,*) "eps2rt_om1",eps2rt_om1," eps2rt_om2",eps2rt_om2,
1592 ! & " eps2rt_om12",eps2rt_om12
1593 ! Calculate whole angle-dependent part of epsilon and contributions
1594 ! to its derivatives
1596 end subroutine sc_angular
1597 !-----------------------------------------------------------------------------
1599 !-----------------------------------------------------------------------------
1600 subroutine int_bounds(total_ints,lower_bound,upper_bound)
1601 ! implicit real*8 (a-h,o-z)
1602 ! include 'DIMENSIONS'
1604 ! include 'COMMON.SETUP'
1605 integer :: total_ints,lower_bound,upper_bound,nint
1606 integer,dimension(0:nfgtasks) :: int4proc,sint4proc !(0:max_fg_procs)
1607 integer :: i,nexcess
1608 nint=total_ints/nfgtasks
1612 nexcess=total_ints-nint*nfgtasks
1614 int4proc(nfgtasks-i)=int4proc(nfgtasks-i)+1
1618 lower_bound=lower_bound+int4proc(i)
1620 upper_bound=lower_bound+int4proc(fg_rank)
1621 lower_bound=lower_bound+1
1623 end subroutine int_bounds
1624 !-----------------------------------------------------------------------------
1625 subroutine int_bounds1(total_ints,lower_bound,upper_bound)
1626 ! implicit real*8 (a-h,o-z)
1627 ! include 'DIMENSIONS'
1629 ! include 'COMMON.SETUP'
1630 integer :: total_ints,lower_bound,upper_bound,nint
1631 integer :: nexcess,i
1632 integer,dimension(0:nfgtasks) :: int4proc,sint4proc !(0:max_fg_procs)
1633 nint=total_ints/nfgtasks1
1637 nexcess=total_ints-nint*nfgtasks1
1639 int4proc(nfgtasks1-i)=int4proc(nfgtasks1-i)+1
1643 lower_bound=lower_bound+int4proc(i)
1645 upper_bound=lower_bound+int4proc(fg_rank1)
1646 lower_bound=lower_bound+1
1648 end subroutine int_bounds1
1649 !-----------------------------------------------------------------------------
1651 !-----------------------------------------------------------------------------
1652 subroutine chainbuild_cart
1653 ! implicit real*8 (a-h,o-z)
1654 ! include 'DIMENSIONS'
1659 ! include 'COMMON.SETUP'
1660 ! include 'COMMON.CHAIN'
1661 ! include 'COMMON.LOCAL'
1662 ! include 'COMMON.TIME1'
1663 ! include 'COMMON.IOUNITS'
1664 integer :: j,i,ierror,ierr
1665 real(kind=8) :: time00,time01
1667 if (nfgtasks.gt.1) then
1668 ! write (iout,*) "BCAST in chainbuild_cart"
1670 ! Broadcast the order to build the chain and compute internal coordinates
1671 ! to the slaves. The slaves receive the order in ERGASTULUM.
1673 ! write (iout,*) "CHAINBUILD_CART: DC before BCAST"
1675 ! write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1676 ! & (dc(j,i+nres),j=1,3)
1679 call MPI_Bcast(7,1,MPI_INTEGER,king,FG_COMM,IERROR)
1680 time_bcast7=time_bcast7+MPI_Wtime()-time00
1682 call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,&
1684 ! write (iout,*) "CHAINBUILD_CART: DC after BCAST"
1686 ! write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1687 ! & (dc(j,i+nres),j=1,3)
1689 ! write (iout,*) "End BCAST in chainbuild_cart"
1691 time_bcast=time_bcast+MPI_Wtime()-time00
1692 time_bcastc=time_bcastc+MPI_Wtime()-time01
1700 c(j,i)=c(j,i-1)+dc(j,i-1)
1705 c(j,i+nres)=c(j,i)+dc(j,i+nres)
1708 ! write (iout,*) "CHAINBUILD_CART"
1710 call int_from_cart1(.false.)
1712 end subroutine chainbuild_cart
1713 !-----------------------------------------------------------------------------
1715 !-----------------------------------------------------------------------------
1716 real(kind=8) function alpha(i1,i2,i3)
1718 ! Calculates the planar angle between atoms (i1), (i2), and (i3).
1720 ! implicit real*8 (a-h,o-z)
1721 ! include 'DIMENSIONS'
1722 ! include 'COMMON.GEO'
1723 ! include 'COMMON.CHAIN'
1726 real(kind=8) :: x12,x23,y12,y23,z12,z23,vnorm,wnorm,scalar
1733 vnorm=dsqrt(x12*x12+y12*y12+z12*z12)
1734 wnorm=dsqrt(x23*x23+y23*y23+z23*z23)
1735 scalar=(x12*x23+y12*y23+z12*z23)/(vnorm*wnorm)
1739 !-----------------------------------------------------------------------------
1740 real(kind=8) function beta(i1,i2,i3,i4)
1742 ! Calculates the dihedral angle between atoms (i1), (i2), (i3) and (i4)
1744 ! implicit real*8 (a-h,o-z)
1745 ! include 'DIMENSIONS'
1746 ! include 'COMMON.GEO'
1747 ! include 'COMMON.CHAIN'
1749 integer :: i1,i2,i3,i4
1750 real(kind=8) :: x12,x23,x34,y12,y23,y34,z12,z23,z34
1751 real(kind=8) :: wx,wy,wz,wnorm,vx,vy,vz,vnorm,scalar,angle
1752 real(kind=8) :: tx,ty,tz
1762 !d print '(2i3,3f10.5)',i1,i2,x12,y12,z12
1763 !d print '(2i3,3f10.5)',i2,i3,x23,y23,z23
1764 !d print '(2i3,3f10.5)',i3,i4,x34,y34,z34
1768 wnorm=dsqrt(wx*wx+wy*wy+wz*wz)
1772 vnorm=dsqrt(vx*vx+vy*vy+vz*vz)
1773 if (vnorm.gt.1.0D-13 .and. wnorm.gt.1.0D-13) then
1774 scalar=(vx*wx+vy*wy+vz*wz)/(vnorm*wnorm)
1775 if (dabs(scalar).gt.1.0D0) &
1776 scalar=0.99999999999999D0*scalar/dabs(scalar)
1778 !d print '(2i4,10f7.3)',i2,i3,vx,vy,vz,wx,wy,wz,vnorm,wnorm,
1783 ! if (angle.le.0.0D0) angle=pi+angle
1787 scalar=tx*x23+ty*y23+tz*z23
1788 if (scalar.lt.0.0D0) angle=-angle
1792 !-----------------------------------------------------------------------------
1793 real(kind=8) function dist(i1,i2)
1795 ! Calculates the distance between atoms (i1) and (i2).
1797 ! implicit real*8 (a-h,o-z)
1798 ! include 'DIMENSIONS'
1799 ! include 'COMMON.GEO'
1800 ! include 'COMMON.CHAIN'
1803 real(kind=8) :: x12,y12,z12
1807 dist=dsqrt(x12*x12+y12*y12+z12*z12)
1811 !-----------------------------------------------------------------------------
1813 !-----------------------------------------------------------------------------
1814 subroutine local_move_init(debug)
1818 ! implicit real*8 (a-h,o-z)
1819 ! include 'DIMENSIONS' ! Needed by COMMON.LOCAL
1820 ! include 'COMMON.GEO' ! For pi, deg2rad
1821 ! include 'COMMON.LOCAL' ! For vbl
1822 ! include 'COMMON.LOCMOVE'
1828 ! Determine wheter to do some debugging output
1829 locmove_output=debug
1831 ! Set the init_called flag to 1
1834 ! The following are never changed
1835 min_theta=60.D0*deg2rad ! (0,PI)
1836 max_theta=175.D0*deg2rad ! (0,PI)
1837 dmin2=vbl*vbl*2.*(1.-cos(min_theta))
1838 dmax2=vbl*vbl*2.*(1.-cos(max_theta))
1841 small2=0.5*small*small
1843 ! Not really necessary...
1849 end subroutine local_move_init
1850 !-----------------------------------------------------------------------------
1851 subroutine local_move(n_start, n_end, PHImin, PHImax)
1852 ! Perform a local move between residues m and n (inclusive)
1853 ! PHImin and PHImax [0,PI] determine the size of the move
1854 ! Works on whatever structure is in the variables theta and phi,
1855 ! sidechain variables are left untouched
1856 ! The final structure is NOT minimized, but both the cartesian
1857 ! variables c and the angles are up-to-date at the end (no further
1858 ! chainbuild is required)
1860 use random,only:ran_number
1862 ! implicit real*8 (a-h,o-z)
1863 ! include 'DIMENSIONS'
1864 ! include 'COMMON.GEO'
1865 ! include 'COMMON.CHAIN'
1866 ! include 'COMMON.VAR'
1867 ! include 'COMMON.MINIM'
1868 ! include 'COMMON.SBRIDGE'
1869 ! include 'COMMON.LOCMOVE'
1871 ! External functions
1872 !EL integer move_res
1873 !EL external move_res
1874 !EL double precision ran_number
1875 !EL external ran_number
1878 integer :: n_start, n_end ! First and last residues to move
1879 real(kind=8) :: PHImin, PHImax ! min/max angles [0,PI]
1883 real(kind=8) :: min,max
1887 ! Check if local_move_init was called. This assumes that it
1888 ! would not be 1 if not explicitely initialized
1889 if (init_called.ne.1) then
1890 write(6,*)' *** local_move_init not called!!!'
1894 ! Quick check for crazy range
1895 if (n_start.gt.n_end .or. n_start.lt.1 .or. n_end.gt.nres) then
1896 write(6,'(a,i3,a,i3)') &
1897 ' *** Cannot make local move between n_start = ',&
1898 n_start,' and n_end = ',n_end
1902 ! Take care of end residues first...
1903 if (n_start.eq.1) then
1904 ! Move residue 1 (completely random)
1905 theta(3)=ran_number(min_theta,max_theta)
1906 phi(4)=ran_number(-PI,PI)
1911 if (n_end.eq.nres) then
1912 ! Move residue nres (completely random)
1913 theta(nres)=ran_number(min_theta,max_theta)
1914 phi(nres)=ran_number(-PI,PI)
1920 ! ...then go through all other residues one by one
1921 ! Start from the two extremes and converge
1926 !$$$c Move the first two residues by less than the others
1927 !$$$ if (i-n_start.lt.3) then
1928 !$$$ if (i-n_start.eq.0) then
1931 !$$$ else if (i-n_start.eq.1) then
1934 !$$$ else if (i-n_start.eq.2) then
1940 ! The actual move, on residue i
1941 iretcode=move_res(min,max,i) ! Discard iretcode
1947 !$$$c Move the last two residues by less than the others
1948 !$$$ if (n_end-j.lt.3) then
1949 !$$$ if (n_end-j.eq.0) then
1952 !$$$ else if (n_end-j.eq.1) then
1955 !$$$ else if (n_end-j.eq.2) then
1961 ! The actual move, on residue j
1962 iretcode=move_res(min,max,j) ! Discard iretcode
1967 call int_from_cart(.false.,.false.)
1970 end subroutine local_move
1971 !-----------------------------------------------------------------------------
1972 subroutine output_tabs
1973 ! Prints out the contents of a_..., b_..., res_...
1977 ! include 'COMMON.GEO'
1978 ! include 'COMMON.LOCMOVE'
1984 write(6,'(8f7.1)')(a_ang(i)*rad2deg,i=0,a_n-1)
1985 write(6,'(8(2x,3l1,2x))')((a_tab(i,j),i=0,2),j=0,a_n-1)
1988 write(6,'(4f7.1)')(b_ang(i)*rad2deg,i=0,b_n-1)
1989 write(6,'(4(2x,3l1,2x))')((b_tab(i,j),i=0,2),j=0,b_n-1)
1992 write(6,'(12f7.1)')(res_ang(i)*rad2deg,i=0,res_n-1)
1993 write(6,'(12(2x,3l1,2x))')((res_tab(0,i,j),i=0,2),j=0,res_n-1)
1994 write(6,'(12(2x,3l1,2x))')((res_tab(1,i,j),i=0,2),j=0,res_n-1)
1995 write(6,'(12(2x,3l1,2x))')((res_tab(2,i,j),i=0,2),j=0,res_n-1)
1998 end subroutine output_tabs
1999 !-----------------------------------------------------------------------------
2000 subroutine angles2tab(PHImin,PHImax,n,ang,tab)
2001 ! Only uses angles if [0,PI] (but PHImin cannot be 0.,
2002 ! and PHImax cannot be PI)
2006 ! include 'COMMON.GEO'
2009 real(kind=8) :: PHImin,PHImax
2013 real(kind=8),dimension(0:3) :: ang
2014 logical,dimension(0:2,0:3) :: tab
2017 if (PHImin .eq. PHImax) then
2018 ! Special case with two 010's
2028 else if (PHImin .eq. PI) then
2029 ! Special case with one 010
2035 else if (PHImax .eq. 0.) then
2036 ! Special case with one 010
2045 if (PHImin .gt. 0.) then
2046 ! Start of range (011)
2051 ! End of range (110)
2055 tab(2,n+1) = .false.
2058 if (PHImax .lt. PI) then
2059 ! Start of range (011)
2064 ! End of range (110)
2068 tab(2,n+1) = .false.
2074 end subroutine angles2tab
2075 !-----------------------------------------------------------------------------
2076 subroutine minmax_angles(x,y,z,r,n,ang,tab)
2077 ! When solutions do not exist, assume all angles
2078 ! are acceptable - i.e., initial geometry must be correct
2082 ! include 'COMMON.GEO'
2083 ! include 'COMMON.LOCMOVE'
2086 real(kind=8) :: x,y,z,r
2090 real(kind=8),dimension(0:3) :: ang
2091 logical,dimension(0:2,0:3) :: tab
2094 real(kind=8) :: num, denom, phi
2095 real(kind=8) :: Kmin, Kmax
2099 num = x*x + y*y + z*z
2102 if (denom .gt. 0.) then
2104 denom = 2.*r*sqrt(denom)
2106 Kmin = (num - dmin2)/denom
2107 Kmax = (num - dmax2)/denom
2109 ! Allowed values of K (else all angles are acceptable)
2112 if (Kmin .gt. 1. .or. abs(Kmin-1.) .lt. small2) then
2114 else if (Kmin .lt. -1. .or. abs(Kmin+1.) .lt. small2) then
2120 if (Kmax .lt. -1. .or. abs(Kmax+1.) .lt. small2) then
2122 else if (Kmax .gt. 1. .or. abs(Kmax-1.) .lt. small2) then
2128 if (Kmax .lt. Kmin) Kmax = Kmin
2130 call angles2tab(Kmin, Kmax, n, ang, tab)
2132 ! Add phi and check that angles are within range (-PI,PI]
2135 if (ang(i) .le. -PI) then
2136 ang(i) = ang(i)+2.*PI
2137 else if (ang(i) .gt. PI) then
2138 ang(i) = ang(i)-2.*PI
2144 end subroutine minmax_angles
2145 !-----------------------------------------------------------------------------
2146 subroutine construct_tab
2147 ! Take a_... and b_... values and produces the results res_...
2148 ! x_ang are assumed to be all different (diff > small)
2149 ! x_tab(1,i) must be 1 for all i (i.e., all x_ang are acceptable)
2153 ! include 'COMMON.LOCMOVE'
2156 integer :: n_max,i,j,index
2162 if (n_max .eq. 0) then
2169 res_tab(j,0,i) = .true.
2170 res_tab(j,2,i) = .true.
2171 res_tab(j,1,i) = .false.
2178 do while (.not.done)
2179 res_ang(index) = flag
2183 if ((a_ang(i)-phi).gt.small .and. &
2184 a_ang(i) .lt. res_ang(index)) then
2185 ! Found a lower angle
2186 res_ang(index) = a_ang(i)
2187 ! Copy the values from a_tab into res_tab(0,,)
2188 res_tab(0,0,index) = a_tab(0,i)
2189 res_tab(0,1,index) = a_tab(1,i)
2190 res_tab(0,2,index) = a_tab(2,i)
2191 ! Set default values for res_tab(1,,)
2192 res_tab(1,0,index) = .true.
2193 res_tab(1,1,index) = .false.
2194 res_tab(1,2,index) = .true.
2195 else if (abs(a_ang(i)-res_ang(index)).lt.small) then
2196 ! Found an equal angle (can only be equal to a b_ang)
2197 res_tab(0,0,index) = a_tab(0,i)
2198 res_tab(0,1,index) = a_tab(1,i)
2199 res_tab(0,2,index) = a_tab(2,i)
2204 if ((b_ang(i)-phi).gt.small .and. &
2205 b_ang(i) .lt. res_ang(index)) then
2206 ! Found a lower angle
2207 res_ang(index) = b_ang(i)
2208 ! Copy the values from b_tab into res_tab(1,,)
2209 res_tab(1,0,index) = b_tab(0,i)
2210 res_tab(1,1,index) = b_tab(1,i)
2211 res_tab(1,2,index) = b_tab(2,i)
2212 ! Set default values for res_tab(0,,)
2213 res_tab(0,0,index) = .true.
2214 res_tab(0,1,index) = .false.
2215 res_tab(0,2,index) = .true.
2216 else if (abs(b_ang(i)-res_ang(index)).lt.small) then
2217 ! Found an equal angle (can only be equal to an a_ang)
2218 res_tab(1,0,index) = b_tab(0,i)
2219 res_tab(1,1,index) = b_tab(1,i)
2220 res_tab(1,2,index) = b_tab(2,i)
2224 if (res_ang(index) .eq. flag) then
2227 else if (index .eq. n_max-1) then
2231 phi = res_ang(index) ! Store previous angle
2239 if (a_n .gt. 0) then
2240 do while (.not.res_tab(0,1,index))
2243 done = res_tab(0,2,index)
2244 do i=index+1,res_n-1
2245 if (res_tab(0,1,i)) then
2246 done = res_tab(0,2,i)
2248 res_tab(0,0,i) = done
2249 res_tab(0,1,i) = done
2250 res_tab(0,2,i) = done
2253 done = res_tab(0,0,index)
2255 if (res_tab(0,1,i)) then
2256 done = res_tab(0,0,i)
2258 res_tab(0,0,i) = done
2259 res_tab(0,1,i) = done
2260 res_tab(0,2,i) = done
2265 res_tab(0,0,i) = .true.
2266 res_tab(0,1,i) = .true.
2267 res_tab(0,2,i) = .true.
2272 if (b_n .gt. 0) then
2273 do while (.not.res_tab(1,1,index))
2276 done = res_tab(1,2,index)
2277 do i=index+1,res_n-1
2278 if (res_tab(1,1,i)) then
2279 done = res_tab(1,2,i)
2281 res_tab(1,0,i) = done
2282 res_tab(1,1,i) = done
2283 res_tab(1,2,i) = done
2286 done = res_tab(1,0,index)
2288 if (res_tab(1,1,i)) then
2289 done = res_tab(1,0,i)
2291 res_tab(1,0,i) = done
2292 res_tab(1,1,i) = done
2293 res_tab(1,2,i) = done
2298 res_tab(1,0,i) = .true.
2299 res_tab(1,1,i) = .true.
2300 res_tab(1,2,i) = .true.
2304 ! Finally fill the last row with AND operation
2307 res_tab(2,j,i) = (res_tab(0,j,i) .and. res_tab(1,j,i))
2312 end subroutine construct_tab
2313 !-----------------------------------------------------------------------------
2314 subroutine construct_ranges(phi_n,phi_start,phi_end)
2315 ! Given the data in res_..., construct a table of
2316 ! min/max allowed angles
2320 ! include 'COMMON.GEO'
2321 ! include 'COMMON.LOCMOVE'
2325 real(kind=8),dimension(0:11) :: phi_start,phi_end
2332 if (res_n .eq. 0) then
2333 ! Any move is allowed
2341 do while (.not.done)
2342 ! Find start of range (01x)
2344 do while (.not.done)
2345 if (res_tab(2,0,index).or.(.not.res_tab(2,1,index))) then
2349 phi_start(phi_n) = res_ang(index)
2351 if (index .eq. res_n) done = .true.
2353 ! If a start was found (index < res_n), find the end of range (x10)
2354 ! It may not be found without wrapping around
2355 if (index .lt. res_n) then
2357 do while (.not.done)
2358 if ((.not.res_tab(2,1,index)).or.res_tab(2,2,index)) then
2363 if (index .eq. res_n) done = .true.
2365 if (index .lt. res_n) then
2366 ! Found the end of the range
2367 phi_end(phi_n) = res_ang(index)
2370 if (index .eq. res_n) then
2376 ! Need to wrap around
2378 phi_end(phi_n) = flag
2382 ! Take care of the last one if need to wrap around
2383 if (phi_end(phi_n) .eq. flag) then
2385 do while ((.not.res_tab(2,1,index)).or.res_tab(2,2,index))
2388 phi_end(phi_n) = res_ang(index) + 2.*PI
2394 end subroutine construct_ranges
2395 !-----------------------------------------------------------------------------
2396 subroutine fix_no_moves(phi)
2400 ! include 'COMMON.GEO'
2401 ! include 'COMMON.LOCMOVE'
2408 real(kind=8) :: diff,temp
2411 ! Look for first 01x in gammas (there MUST be at least one)
2414 do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index)))
2417 if (res_ang(index) .le. 0.D0) then ! Make sure it's from PHImax
2418 ! Try to increase PHImax
2419 if (index .gt. 0) then
2420 phi = res_ang(index-1)
2421 diff = abs(res_ang(index) - res_ang(index-1))
2423 ! Look for last (corresponding) x10
2425 do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index))
2428 if (index .lt. res_n-1) then
2429 temp = abs(res_ang(index) - res_ang(index+1))
2430 if (temp .lt. diff) then
2431 phi = res_ang(index+1)
2437 ! If increasing PHImax didn't work, decreasing PHImin
2438 ! will (with one exception)
2439 ! Look for first x10 (there MUST be at least one)
2441 do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index))
2444 if (res_ang(index) .lt. 0.D0) then ! Make sure it's from PHImin
2445 ! Try to decrease PHImin
2446 if (index .lt. res_n-1) then
2447 temp = abs(res_ang(index) - res_ang(index+1))
2448 if (res_ang(index+1) .le. 0.D0 .and. temp .lt. diff) then
2449 phi = res_ang(index+1)
2453 ! Look for last (corresponding) 01x
2455 do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index)))
2458 if (index .gt. 0) then
2459 temp = abs(res_ang(index) - res_ang(index-1))
2460 if (res_ang(index-1) .ge. 0.D0 .and. temp .lt. diff) then
2461 phi = res_ang(index-1)
2467 ! If it still didn't work, it must be PHImax == 0. or PHImin == PI
2468 if (diff .eq. flag) then
2470 if (res_tab(index,1,0) .or. (.not.res_tab(index,1,1)) .or. &
2471 res_tab(index,1,2)) index = res_n - 1
2472 ! This MUST work at this point
2473 if (index .eq. 0) then
2476 phi = res_ang(index - 1)
2481 end subroutine fix_no_moves
2482 !-----------------------------------------------------------------------------
2483 integer function move_res(PHImin,PHImax,i_move)
2484 ! Moves residue i_move (in array c), leaving everything else fixed
2485 ! Starting geometry is not checked, it should be correct!
2486 ! R(,i_move) is the only residue that will move, but must have
2487 ! 1 < i_move < nres (i.e., cannot move ends)
2488 ! Whether any output is done is controlled by locmove_output
2490 use random,only:ran_number
2492 ! implicit real*8 (a-h,o-z)
2493 ! include 'DIMENSIONS'
2494 ! include 'COMMON.CHAIN'
2495 ! include 'COMMON.GEO'
2496 ! include 'COMMON.LOCMOVE'
2498 ! External functions
2499 !EL double precision ran_number
2500 !EL external ran_number
2503 real(kind=8) :: PHImin,PHImax
2507 ! 0: move successfull
2508 ! 1: Dmin or Dmax had to be modified
2509 ! 2: move failed - check your input geometry
2513 real(kind=8),dimension(0:2) :: X,Y,Z,Orig
2514 real(kind=8),dimension(0:2) :: P
2515 logical :: no_moves,done
2516 integer :: index,i,j
2517 real(kind=8) :: phi,temp,radius
2518 real(kind=8),dimension(0:11) :: phi_start,phi_end
2521 ! Set up the coordinate system
2523 Orig(i)=0.5*(c(i+1,i_move-1)+c(i+1,i_move+1)) ! Position of origin
2527 Z(i)=c(i+1,i_move+1)-c(i+1,i_move-1)
2529 temp=sqrt(Z(0)*Z(0)+Z(1)*Z(1)+Z(2)*Z(2))
2535 X(i)=c(i+1,i_move)-Orig(i)
2537 ! radius is the radius of the circle on which c(,i_move) can move
2538 radius=sqrt(X(0)*X(0)+X(1)*X(1)+X(2)*X(2))
2543 Y(0)=Z(1)*X(2)-X(1)*Z(2)
2544 Y(1)=X(0)*Z(2)-Z(0)*X(2)
2545 Y(2)=Z(0)*X(1)-X(0)*Z(1)
2547 ! Calculate min, max angles coming from dmin, dmax to c(,i_move-2)
2548 if (i_move.gt.2) then
2550 P(i)=c(i+1,i_move-2)-Orig(i)
2552 call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2),&
2553 P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2),&
2554 P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2),&
2555 radius,a_n,a_ang,a_tab)
2560 ! Calculate min, max angles coming from dmin, dmax to c(,i_move+2)
2561 if (i_move.lt.nres-2) then
2563 P(i)=c(i+1,i_move+2)-Orig(i)
2565 call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2),&
2566 P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2),&
2567 P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2),&
2568 radius,b_n,b_ang,b_tab)
2573 ! Construct the resulting table for alpha and beta
2574 call construct_tab()
2576 if (locmove_output) then
2577 print *,'ALPHAS & BETAS TABLE'
2581 ! Check that there is at least one possible move
2583 if (res_n .eq. 0) then
2587 do while ((index .lt. res_n) .and. no_moves)
2588 if (res_tab(2,1,index)) no_moves = .false.
2593 if (locmove_output) print *,' *** Cannot move anywhere'
2598 ! Transfer res_... into a_...
2601 if ( (res_tab(2,0,i).neqv.res_tab(2,1,i)) .or. &
2602 (res_tab(2,0,i).neqv.res_tab(2,2,i)) ) then
2603 a_ang(a_n) = res_ang(i)
2605 a_tab(j,a_n) = res_tab(2,j,i)
2611 ! Check that the PHI's are within [0,PI]
2612 if (PHImin .lt. 0. .or. abs(PHImin) .lt. small) PHImin = -flag
2613 if (PHImin .gt. PI .or. abs(PHImin-PI) .lt. small) PHImin = PI
2614 if (PHImax .gt. PI .or. abs(PHImax-PI) .lt. small) PHImax = flag
2615 if (PHImax .lt. 0. .or. abs(PHImax) .lt. small) PHImax = 0.
2616 if (PHImax .lt. PHImin) PHImax = PHImin
2617 ! Calculate min and max angles coming from PHImin and PHImax,
2618 ! and put them in b_...
2619 call angles2tab(PHImin, PHImax, b_n, b_ang, b_tab)
2620 ! Construct the final table
2621 call construct_tab()
2623 if (locmove_output) then
2624 print *,'FINAL TABLE'
2628 ! Check that there is at least one possible move
2630 if (res_n .eq. 0) then
2634 do while ((index .lt. res_n) .and. no_moves)
2635 if (res_tab(2,1,index)) no_moves = .false.
2641 ! Take care of the case where no solution exists...
2642 call fix_no_moves(phi)
2643 if (locmove_output) then
2644 print *,' *** Had to modify PHImin or PHImax'
2645 print *,'phi: ',phi*rad2deg
2649 ! ...or calculate the solution
2650 ! Construct phi_start/phi_end arrays
2651 call construct_ranges(phi_n, phi_start, phi_end)
2652 ! Choose random angle phi in allowed range(s)
2655 temp = temp + phi_end(i) - phi_start(i)
2657 phi = ran_number(phi_start(0),phi_start(0)+temp)
2660 do while (.not.done)
2661 if (phi .lt. phi_end(index)) then
2666 if (index .eq. phi_n) then
2668 else if (.not.done) then
2669 phi = phi + phi_start(index) - phi_end(index-1)
2672 if (index.eq.phi_n) phi=phi_end(phi_n-1) ! Fix numerical errors
2673 if (phi .gt. PI) phi = phi-2.*PI
2675 if (locmove_output) then
2676 print *,'ALLOWED RANGE(S)'
2678 print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg
2680 print *,'phi: ',phi*rad2deg
2685 ! Re-use radius as temp variable
2686 temp=radius*cos(phi)
2687 radius=radius*sin(phi)
2689 c(i+1,i_move)=Orig(i)+temp*X(i)+radius*Y(i)
2693 end function move_res
2694 !-----------------------------------------------------------------------------
2699 ! implicit real*8 (a-h,o-z)
2700 ! include 'DIMENSIONS'
2701 ! include 'COMMON.GEO'
2702 ! include 'COMMON.LOCAL'
2703 ! include 'COMMON.LOCMOVE'
2705 ! External functions
2706 !EL integer move_res
2707 !EL external move_res
2712 real(kind=8),dimension(0:11) :: phi_start,phi_end
2714 real(kind=8),dimension(0:2,0:5) :: R
2716 locmove_output=.true.
2718 ! call angles2tab(30.*deg2rad,70.*deg2rad,a_n,a_ang,a_tab)
2719 ! call angles2tab(80.*deg2rad,130.*deg2rad,b_n,b_ang,b_tab)
2720 ! call minmax_angles(0.D0,3.8D0,0.D0,3.8D0,b_n,b_ang,b_tab)
2721 ! call construct_tab
2724 ! call construct_ranges(phi_n,phi_start,phi_end)
2726 ! print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg
2729 ! call fix_no_moves(phi)
2730 ! print *,'NO MOVES FOUND, BEST PHI IS',phi*rad2deg
2736 R(1,1)=-cos(28.D0*deg2rad)
2737 R(2,1)=-0.5D0-sin(28.D0*deg2rad)
2741 R(0,3)=cos(30.D0*deg2rad)
2748 R(1,5)=cos(26.D0*deg2rad)
2749 R(2,5)=0.5D0+sin(26.D0*deg2rad)
2755 ! i=move_res(R(0,1),0.D0*deg2rad,180.D0*deg2rad)
2757 i=move_res(0.D0*deg2rad,180.D0*deg2rad,imov)
2758 print *,'RETURNED ',i
2759 print *,(R(i,3)/vbl,i=0,2)
2762 end subroutine loc_test
2764 !-----------------------------------------------------------------------------
2766 !-----------------------------------------------------------------------------
2767 subroutine MATMULT(A1,A2,A3)
2768 ! implicit real*8 (a-h,o-z)
2769 ! include 'DIMENSIONS'
2772 real(kind=8) :: A3IJ
2774 real(kind=8),DIMENSION(3,3) :: A1,A2,A3
2775 real(kind=8),DIMENSION(3,3) :: AI3
2780 3 A3IJ=A3IJ+A1(I,K)*A2(K,J)
2788 end subroutine MATMULT
2789 !-----------------------------------------------------------------------------
2791 !-----------------------------------------------------------------------------
2792 subroutine int_from_cart(lside,lprn)
2793 ! implicit real*8 (a-h,o-z)
2794 ! include 'DIMENSIONS'
2795 use control_data,only:out1file
2799 ! include 'COMMON.LOCAL'
2800 ! include 'COMMON.VAR'
2801 ! include 'COMMON.CHAIN'
2802 ! include 'COMMON.INTERACT'
2803 ! include 'COMMON.IOUNITS'
2804 ! include 'COMMON.GEO'
2805 ! include 'COMMON.NAMES'
2806 ! include 'COMMON.CONTROL'
2807 ! include 'COMMON.SETUP'
2808 character(len=3) :: seq,res
2810 character(len=80) :: card
2811 real(kind=8),dimension(3,20) :: sccor
2812 integer :: i,j,iti !el rescode,
2813 logical :: lside,lprn
2814 real(kind=8) :: di,cosfac,sinfac
2818 if(me.eq.king.or..not.out1file)then
2820 write (iout,'(/a)') &
2821 'Internal coordinates calculated from crystal structure.'
2823 write (iout,'(8a)') ' Res ',' dvb',' Theta',&
2824 ' Gamma',' Dsc_id',' Dsc',' Alpha',&
2827 write (iout,'(4a)') ' Res ',' dvb',' Theta',&
2834 if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
2835 write (iout,'(a,i4)') 'Bad Cartesians for residue',i
2838 vbld(i+1)=dist(i,i+1)
2839 vbld_inv(i+1)=1.0d0/vbld(i+1)
2840 if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
2841 if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
2843 ! if (unres_pdb) then
2844 ! if (itype(1).eq.21) then
2845 ! theta(3)=90.0d0*deg2rad
2846 ! phi(4)=180.0d0*deg2rad
2848 ! vbld_inv(2)=1.0d0/vbld(2)
2850 ! if (itype(nres).eq.21) then
2851 ! theta(nres)=90.0d0*deg2rad
2852 ! phi(nres)=180.0d0*deg2rad
2854 ! vbld_inv(nres)=1.0d0/vbld(2)
2860 c(j,nres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i) &
2861 +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
2865 ! 10/03/12 Adam: Correction for zero SC-SC bond length
2866 if (itype(i).ne.10 .and. itype(i).ne.ntyp1 .and. di.eq.0.0d0) &
2869 if (itype(i).ne.10) then
2870 vbld_inv(i+nres)=1.0d0/di
2872 vbld_inv(i+nres)=0.0d0
2875 alph(i)=alpha(nres+i,i,nres2)
2876 omeg(i)=beta(nres+i,i,nres2,i+1)
2878 if(me.eq.king.or..not.out1file)then
2880 write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),&
2881 rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),&
2882 rad2deg*alph(i),rad2deg*omeg(i)
2888 if(me.eq.king.or..not.out1file) &
2889 write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),&
2890 rad2deg*theta(i),rad2deg*phi(i)
2894 end subroutine int_from_cart
2895 !-----------------------------------------------------------------------------
2896 subroutine sc_loc_geom(lprn)
2897 ! implicit real*8 (a-h,o-z)
2898 ! include 'DIMENSIONS'
2899 use control_data,only:out1file
2903 ! include 'COMMON.LOCAL'
2904 ! include 'COMMON.VAR'
2905 ! include 'COMMON.CHAIN'
2906 ! include 'COMMON.INTERACT'
2907 ! include 'COMMON.IOUNITS'
2908 ! include 'COMMON.GEO'
2909 ! include 'COMMON.NAMES'
2910 ! include 'COMMON.CONTROL'
2911 ! include 'COMMON.SETUP'
2912 real(kind=8),dimension(3) :: x_prime,y_prime,z_prime
2915 integer :: i,j,it,iti
2916 real(kind=8) :: cosfac2,sinfac2,xx,yy,zz,cosfac,sinfac
2919 dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
2923 if (itype(i).ne.10) then
2925 dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
2929 dc_norm(j,i+nres)=0.0d0
2934 costtab(i+1) =dcos(theta(i+1))
2935 sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
2936 cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
2937 sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
2938 cosfac2=0.5d0/(1.0d0+costtab(i+1))
2939 cosfac=dsqrt(cosfac2)
2940 sinfac2=0.5d0/(1.0d0-costtab(i+1))
2941 sinfac=dsqrt(sinfac2)
2945 ! Compute the axes of tghe local cartesian coordinates system; store in
2946 ! x_prime, y_prime and z_prime
2954 x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
2955 y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
2957 call vecpr(x_prime,y_prime,z_prime)
2959 ! Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
2960 ! to local coordinate system. Store in xx, yy, zz.
2966 xx = xx + x_prime(j)*dc_norm(j,i+nres)
2967 yy = yy + y_prime(j)*dc_norm(j,i+nres)
2968 zz = zz + z_prime(j)*dc_norm(j,i+nres)
2983 if(me.eq.king.or..not.out1file) &
2984 write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),&
2989 end subroutine sc_loc_geom
2990 !-----------------------------------------------------------------------------
2991 subroutine sccenter(ires,nscat,sccor)
2992 ! implicit real*8 (a-h,o-z)
2993 ! include 'DIMENSIONS'
2994 ! include 'COMMON.CHAIN'
2995 integer :: i,j,ires,nscat
2996 real(kind=8),dimension(3,20) :: sccor
2997 real(kind=8) :: sccmj
3001 sccmj=sccmj+sccor(j,i)
3003 dc(j,ires)=sccmj/nscat
3006 end subroutine sccenter
3008 !-----------------------------------------------------------------------------
3009 subroutine bond_regular
3011 ! implicit real*8 (a-h,o-z)
3012 ! include 'DIMENSIONS'
3013 ! include 'COMMON.VAR'
3014 ! include 'COMMON.LOCAL'
3015 ! include 'COMMON.CALC'
3016 ! include 'COMMON.INTERACT'
3017 ! include 'COMMON.CHAIN'
3020 vbld_inv(i+1)=1.0d0/vbld(i+1)
3021 vbld(i+1+nres)=dsc(itype(i+1))
3022 vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
3023 ! print *,vbld(i+1),vbld(i+1+nres)
3026 end subroutine bond_regular
3028 !-----------------------------------------------------------------------------
3030 !-----------------------------------------------------------------------------
3031 subroutine refsys(i2,i3,i4,e1,e2,e3,fail)
3032 ! This subroutine calculates unit vectors of a local reference system
3033 ! defined by atoms (i2), (i3), and (i4). The x axis is the axis from
3034 ! atom (i3) to atom (i2), and the xy plane is the plane defined by atoms
3035 ! (i2), (i3), and (i4). z axis is directed according to the sign of the
3036 ! vector product (i3)-(i2) and (i3)-(i4). Sets fail to .true. if atoms
3037 ! (i2) and (i3) or (i3) and (i4) coincide or atoms (i2), (i3), and (i4)
3038 ! form a linear fragment. Returns vectors e1, e2, and e3.
3039 ! implicit real*8 (a-h,o-z)
3040 ! include 'DIMENSIONS'
3042 real(kind=8),dimension(3) :: e1,e2,e3
3043 real(kind=8),dimension(3) :: u,z
3044 ! include 'COMMON.IOUNITS'
3045 ! include 'COMMON.CHAIN'
3046 real(kind=8) :: coinc=1.0D-13,align=1.0D-13
3048 integer :: i,i1,i2,i3,i4
3049 real(kind=8) :: v1,v2,v3,s1,s2,zi,ui,anorm
3062 if (s1.gt.coinc) goto 2
3063 write (iout,1000) i2,i3,i1
3068 2 if (s2.gt.coinc) goto 4
3069 write(iout,1000) i3,i4,i1
3076 v1=z(2)*u(3)-z(3)*u(2)
3077 v2=z(3)*u(1)-z(1)*u(3)
3078 v3=z(1)*u(2)-z(2)*u(1)
3079 anorm=dsqrt(v1*v1+v2*v2+v3*v3)
3080 if (anorm.gt.align) goto 6
3081 write (iout,1010) i2,i3,i4,i1
3093 e2(1)=e1(3)*e3(2)-e1(2)*e3(3)
3094 e2(2)=e1(1)*e3(3)-e1(3)*e3(1)
3095 e2(3)=e1(2)*e3(1)-e1(1)*e3(2)
3096 1000 format (/1x,' * * * Error - atoms',i4,' and',i4,' coincide.',&
3097 'coordinates of atom',i4,' are set to zero.')
3098 1010 format (/1x,' * * * Error - atoms',2(i4,2h, ),i4,' form a linear',&
3099 ' fragment. coordinates of atom',i4,' are set to zero.')
3101 end subroutine refsys
3102 !-----------------------------------------------------------------------------
3104 !-----------------------------------------------------------------------------
3105 subroutine int_to_cart
3106 !--------------------------------------------------------------
3107 ! This subroutine converts the energy derivatives from internal
3108 ! coordinates to cartesian coordinates
3109 !-------------------------------------------------------------
3110 ! implicit real*8 (a-h,o-z)
3111 ! include 'DIMENSIONS'
3112 ! include 'COMMON.VAR'
3113 ! include 'COMMON.CHAIN'
3114 ! include 'COMMON.DERIV'
3115 ! include 'COMMON.GEO'
3116 ! include 'COMMON.LOCAL'
3117 ! include 'COMMON.INTERACT'
3118 ! include 'COMMON.MD'
3119 ! include 'COMMON.IOUNITS'
3120 ! include 'COMMON.SCCOR'
3121 ! calculating dE/ddc1
3124 if (nres.lt.3) go to 18
3126 gcart(j,1)=gcart(j,1)+gloc(1,icg)*dphi(j,1,4) &
3127 +gloc(nres-2,icg)*dtheta(j,1,3)
3128 if(itype(2).ne.10) then
3129 gcart(j,1)=gcart(j,1)+gloc(ialph(2,1),icg)*dalpha(j,1,2)+ &
3130 gloc(ialph(2,1)+nside,icg)*domega(j,1,2)
3133 ! Calculating the remainder of dE/ddc2
3135 gcart(j,2)=gcart(j,2)+gloc(1,icg)*dphi(j,2,4)+ &
3136 gloc(nres-2,icg)*dtheta(j,2,3)+gloc(nres-1,icg)*dtheta(j,1,4)
3137 if(itype(2).ne.10) then
3138 gcart(j,2)=gcart(j,2)+gloc(ialph(2,1),icg)*dalpha(j,2,2)+ &
3139 gloc(ialph(2,1)+nside,icg)*domega(j,2,2)
3141 if(itype(3).ne.10) then
3142 gcart(j,2)=gcart(j,2)+gloc(ialph(3,1),icg)*dalpha(j,1,3)+ &
3143 gloc(ialph(3,1)+nside,icg)*domega(j,1,3)
3146 gcart(j,2)=gcart(j,2)+gloc(2,icg)*dphi(j,1,5)
3149 ! If there are only five residues
3152 gcart(j,3)=gcart(j,3)+gloc(1,icg)*dphi(j,3,4)+gloc(2,icg)* &
3153 dphi(j,2,5)+gloc(nres-1,icg)*dtheta(j,2,4)+gloc(nres,icg)* &
3155 if(itype(3).ne.10) then
3156 gcart(j,3)=gcart(j,3)+gloc(ialph(3,1),icg)* &
3157 dalpha(j,2,3)+gloc(ialph(3,1)+nside,icg)*domega(j,2,3)
3159 if(itype(4).ne.10) then
3160 gcart(j,3)=gcart(j,3)+gloc(ialph(4,1),icg)* &
3161 dalpha(j,1,4)+gloc(ialph(4,1)+nside,icg)*domega(j,1,4)
3165 ! If there are more than five residues
3169 gcart(j,i)=gcart(j,i)+gloc(i-2,icg)*dphi(j,3,i+1) &
3170 +gloc(i-1,icg)*dphi(j,2,i+2)+ &
3171 gloc(i,icg)*dphi(j,1,i+3)+gloc(nres+i-4,icg)*dtheta(j,2,i+1)+ &
3172 gloc(nres+i-3,icg)*dtheta(j,1,i+2)
3173 if(itype(i).ne.10) then
3174 gcart(j,i)=gcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,2,i)+ &
3175 gloc(ialph(i,1)+nside,icg)*domega(j,2,i)
3177 if(itype(i+1).ne.10) then
3178 gcart(j,i)=gcart(j,i)+gloc(ialph(i+1,1),icg)*dalpha(j,1,i+1) &
3179 +gloc(ialph(i+1,1)+nside,icg)*domega(j,1,i+1)
3184 ! Setting dE/ddnres-2
3187 gcart(j,nres-2)=gcart(j,nres-2)+gloc(nres-4,icg)* &
3188 dphi(j,3,nres-1)+gloc(nres-3,icg)*dphi(j,2,nres) &
3189 +gloc(2*nres-6,icg)* &
3190 dtheta(j,2,nres-1)+gloc(2*nres-5,icg)*dtheta(j,1,nres)
3191 if(itype(nres-2).ne.10) then
3192 gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-2,1),icg)* &
3193 dalpha(j,2,nres-2)+gloc(ialph(nres-2,1)+nside,icg)* &
3196 if(itype(nres-1).ne.10) then
3197 gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-1,1),icg)* &
3198 dalpha(j,1,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
3203 ! Settind dE/ddnres-1
3205 gcart(j,nres-1)=gcart(j,nres-1)+gloc(nres-3,icg)*dphi(j,3,nres)+ &
3206 gloc(2*nres-5,icg)*dtheta(j,2,nres)
3207 if(itype(nres-1).ne.10) then
3208 gcart(j,nres-1)=gcart(j,nres-1)+gloc(ialph(nres-1,1),icg)* &
3209 dalpha(j,2,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
3213 ! The side-chain vector derivatives
3215 if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3217 gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i) &
3218 +gloc(ialph(i,1)+nside,icg)*domega(j,3,i)
3222 !----------------------------------------------------------------------
3223 ! INTERTYP=1 SC...Ca...Ca...Ca
3224 ! INTERTYP=2 Ca...Ca...Ca...SC
3225 ! INTERTYP=3 SC...Ca...Ca...SC
3226 ! calculating dE/ddc1
3230 ! write (iout,*) "poczotkoawy",i,gloc_sc(1,i,icg)
3232 if (nres.lt.2) return
3233 if ((nres.lt.3).and.(itype(1).eq.10)) return
3234 if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then
3236 !c Derviative was calculated for oposite vector of side chain therefore
3237 ! there is "-" sign before gloc_sc
3238 gxcart(j,1)=gxcart(j,1)-gloc_sc(1,0,icg)* &
3240 gcart(j,1)=gcart(j,1)+gloc_sc(1,0,icg)* &
3242 if ((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then
3243 gxcart(j,1)= gxcart(j,1) &
3244 -gloc_sc(3,0,icg)*dtauangle(j,3,1,3)
3245 gcart(j,1)=gcart(j,1)+gloc_sc(3,0,icg)* &
3250 if ((nres.ge.3).and.(itype(3).ne.10).and.(itype(3).ne.ntyp1)) &
3253 gcart(j,1)=gcart(j,1)+gloc_sc(2,1,icg)*dtauangle(j,2,1,4)
3256 ! As potetnial DO NOT depend on omicron anlge their derivative is
3258 ! & +gloc_sc(intertyp,nres-2,icg)*dtheta(j,1,3)
3260 ! Calculating the remainder of dE/ddc2
3262 if((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then
3263 if (itype(1).ne.10) gxcart(j,2)=gxcart(j,2)+ &
3264 gloc_sc(3,0,icg)*dtauangle(j,3,3,3)
3265 if ((itype(3).ne.10).and.(nres.ge.3).and.(itype(3).ne.ntyp1)) &
3267 gxcart(j,2)=gxcart(j,2)-gloc_sc(3,1,icg)*dtauangle(j,3,1,4)
3268 !c the - above is due to different vector direction
3269 gcart(j,2)=gcart(j,2)+gloc_sc(3,1,icg)*dtauangle(j,3,2,4)
3272 gxcart(j,2)=gxcart(j,2)-gloc_sc(1,1,icg)*dtauangle(j,1,1,4)
3273 !c the - above is due to different vector direction
3274 gcart(j,2)=gcart(j,2)+gloc_sc(1,1,icg)*dtauangle(j,1,2,4)
3275 ! write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,2,4),"gcart"
3276 ! write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx"
3279 if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then
3280 gcart(j,2)=gcart(j,2)+gloc_sc(1,0,icg)*dtauangle(j,1,3,3)
3281 ! write(iout,*) gloc_sc(1,0,icg),dtauangle(j,1,3,3)
3283 if ((itype(3).ne.10).and.(nres.ge.3)) then
3284 gcart(j,2)=gcart(j,2)+gloc_sc(2,1,icg)*dtauangle(j,2,2,4)
3285 ! write(iout,*) gloc_sc(2,1,icg),dtauangle(j,2,2,4)
3287 if ((itype(4).ne.10).and.(nres.ge.4)) then
3288 gcart(j,2)=gcart(j,2)+gloc_sc(2,2,icg)*dtauangle(j,2,1,5)
3289 ! write(iout,*) gloc_sc(2,2,icg),dtauangle(j,2,1,5)
3292 ! write(iout,*) gcart(j,2),itype(2),itype(1),itype(3), "gcart2"
3294 ! If there are more than five residues
3298 ! write(iout,*) "before", gcart(j,i)
3299 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
3300 gxcart(j,i)=gxcart(j,i)+gloc_sc(2,i-2,icg) &
3301 *dtauangle(j,2,3,i+1) &
3302 -gloc_sc(1,i-1,icg)*dtauangle(j,1,1,i+2)
3303 gcart(j,i)=gcart(j,i)+gloc_sc(1,i-1,icg) &
3304 *dtauangle(j,1,2,i+2)
3305 ! write(iout,*) "new",j,i,
3306 ! & gcart(j,i),gloc_sc(1,i-1,icg),dtauangle(j,1,2,i+2)
3307 if (itype(i-1).ne.10) then
3308 gxcart(j,i)=gxcart(j,i)+gloc_sc(3,i-2,icg) &
3309 *dtauangle(j,3,3,i+1)
3311 if (itype(i+1).ne.10) then
3312 gxcart(j,i)=gxcart(j,i)-gloc_sc(3,i-1,icg) &
3313 *dtauangle(j,3,1,i+2)
3314 gcart(j,i)=gcart(j,i)+gloc_sc(3,i-1,icg) &
3315 *dtauangle(j,3,2,i+2)
3318 if (itype(i-1).ne.10) then
3319 gcart(j,i)=gcart(j,i)+gloc_sc(1,i-2,icg)* &
3320 dtauangle(j,1,3,i+1)
3322 if (itype(i+1).ne.10) then
3323 gcart(j,i)=gcart(j,i)+gloc_sc(2,i-1,icg)* &
3324 dtauangle(j,2,2,i+2)
3325 ! write(iout,*) "numer",i,gloc_sc(2,i-1,icg),
3326 ! & dtauangle(j,2,2,i+2)
3328 if (itype(i+2).ne.10) then
3329 gcart(j,i)=gcart(j,i)+gloc_sc(2,i,icg)* &
3330 dtauangle(j,2,1,i+3)
3335 ! Setting dE/ddnres-1
3338 if ((itype(nres-1).ne.10).and.(itype(nres-1).ne.ntyp1)) then
3339 gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(2,nres-3,icg) &
3340 *dtauangle(j,2,3,nres)
3341 ! write (iout,*) "gxcart(nres-1)", gloc_sc(2,nres-3,icg),
3342 ! & dtauangle(j,2,3,nres), gxcart(j,nres-1)
3343 if (itype(nres-2).ne.10) then
3344 gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) &
3345 *dtauangle(j,3,3,nres)
3347 if ((itype(nres).ne.10).and.(itype(nres).ne.ntyp1)) then
3348 gxcart(j,nres-1)=gxcart(j,nres-1)-gloc_sc(3,nres-2,icg) &
3349 *dtauangle(j,3,1,nres+1)
3350 gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(3,nres-2,icg) &
3351 *dtauangle(j,3,2,nres+1)
3354 if ((itype(nres-2).ne.10).and.(itype(nres-2).ne.ntyp1)) then
3355 gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(1,nres-3,icg)* &
3356 dtauangle(j,1,3,nres)
3358 if ((itype(nres).ne.10).and.(itype(nres).ne.ntyp1)) then
3359 gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(2,nres-2,icg)* &
3360 dtauangle(j,2,2,nres+1)
3361 ! write (iout,*) "gcart(nres-1)", gloc_sc(2,nres-2,icg),
3362 ! & dtauangle(j,2,2,nres+1), itype(nres-1),itype(nres)
3367 if ((nres.ge.3).and.(itype(nres).ne.10).and. &
3368 (itype(nres).ne.ntyp1))then
3370 gxcart(j,nres)=gxcart(j,nres)+gloc_sc(3,nres-2,icg) &
3371 *dtauangle(j,3,3,nres+1)+gloc_sc(2,nres-2,icg) &
3372 *dtauangle(j,2,3,nres+1)
3375 ! The side-chain vector derivatives
3377 end subroutine int_to_cart
3379 !-----------------------------------------------------------------------------
3381 !-----------------------------------------------------------------------------
3382 subroutine gen_dist_constr
3383 ! Generate CA distance constraints.
3384 ! implicit real*8 (a-h,o-z)
3385 ! include 'DIMENSIONS'
3386 ! include 'COMMON.IOUNITS'
3387 ! include 'COMMON.GEO'
3388 ! include 'COMMON.VAR'
3389 ! include 'COMMON.INTERACT'
3390 ! include 'COMMON.LOCAL'
3391 ! include 'COMMON.NAMES'
3392 ! include 'COMMON.CHAIN'
3393 ! include 'COMMON.FFIELD'
3394 ! include 'COMMON.SBRIDGE'
3395 ! include 'COMMON.HEADER'
3396 ! include 'COMMON.CONTROL'
3397 ! include 'COMMON.DBASE'
3398 ! include 'COMMON.THREAD'
3399 ! include 'COMMON.TIME1'
3400 ! integer :: itype_pdb !(maxres)
3401 ! common /pizda/ itype_pdb(nres)
3402 character(len=2) :: iden
3405 !d print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
3406 !d write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
3407 !d & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
3409 do i=nstart_sup,nstart_sup+nsup-1
3410 !d write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
3411 !d & ' seq_pdb', restyp(itype_pdb(i))
3412 do j=i+2,nstart_sup+nsup-1
3414 ihpb(nhpb)=i+nstart_seq-nstart_sup
3415 jhpb(nhpb)=j+nstart_seq-nstart_sup
3417 dhpb(nhpb)=dist(i,j)
3420 !d write (iout,'(a)') 'Distance constraints:'
3425 !d if (ii.gt.nres) then
3430 !d write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
3431 !d & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
3432 !d & dhpb(i),forcon(i)
3434 deallocate(itype_pdb)
3437 end subroutine gen_dist_constr
3439 !-----------------------------------------------------------------------------
3441 !-----------------------------------------------------------------------------
3442 subroutine cartprint
3444 use geometry_data, only: c
3445 use energy_data, only: itype
3446 ! implicit real*8 (a-h,o-z)
3447 ! include 'DIMENSIONS'
3448 ! include 'COMMON.CHAIN'
3449 ! include 'COMMON.INTERACT'
3450 ! include 'COMMON.NAMES'
3451 ! include 'COMMON.IOUNITS'
3456 write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
3457 c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
3459 100 format (//' alpha-carbon coordinates ',&
3460 ' centroid coordinates'/ &
3461 ' ', 6X,'X',11X,'Y',11X,'Z',&
3462 10X,'X',11X,'Y',11X,'Z')
3463 110 format (a,'(',i3,')',6f12.5)
3465 end subroutine cartprint
3466 !-----------------------------------------------------------------------------
3467 !-----------------------------------------------------------------------------
3468 subroutine alloc_geo_arrays
3469 !EL Allocation of tables used by module energy
3471 integer :: i,j,nres2
3475 allocate(phibound(2,nres+2)) !(2,maxres)
3476 !----------------------
3478 ! common /chain/ in molread
3479 ! real(kind=8),dimension(:,:),allocatable :: c !(3,maxres2+2)
3480 ! real(kind=8),dimension(:,:),allocatable :: dc
3481 allocate(dc_old(3,0:nres2))
3482 if(.not.allocated(dc_norm2)) allocate(dc_norm2(3,0:nres2)) !(3,0:maxres2)
3483 !el if(.not.allocated(dc_norm))
3484 !elwrite(iout,*) "jestem w alloc geo 1"
3485 if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:nres2)) !(3,0:maxres2)
3486 !elwrite(iout,*) "jestem w alloc geo 1"
3487 allocate(xloc(3,nres),xrot(3,nres))
3488 !elwrite(iout,*) "jestem w alloc geo 1"
3494 !elwrite(iout,*) "jestem w alloc geo 1"
3495 allocate(dc_work(6*nres)) !(MAXRES6) maxres6=6*maxres
3497 allocate(t(3,3,nres),r(3,3,nres))
3498 allocate(prod(3,3,nres),rt(3,3,nres)) !(3,3,maxres)
3499 ! common /refstruct/
3500 if(.not.allocated(cref)) allocate(cref(3,nres2+2,maxperm)) !(3,maxres2+2,maxperm)
3501 !elwrite(iout,*) "jestem w alloc geo 2"
3502 allocate(crefjlee(3,nres2+2)) !(3,maxres2+2)
3503 if(.not.allocated(chain_rep)) allocate(chain_rep(3,nres2+2,maxsym)) !(3,maxres2+2,maxsym)
3504 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3505 ! common /from_zscore/ in module.compare
3506 !----------------------
3508 ! Inverses of the actual virtual bond lengths
3509 ! common /invlen/ in io_conf: molread or readpdb
3510 ! real(kind=8),dimension(:),allocatable :: vbld_inv !(maxres2)
3511 !----------------------
3513 ! Store the geometric variables in the following COMMON block.
3514 ! common /var/ in readpdb or ...
3515 if(.not.allocated(theta)) allocate(theta(nres+2))
3516 if(.not.allocated(phi)) allocate(phi(nres+2))
3517 if(.not.allocated(alph)) allocate(alph(nres+2))
3518 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3519 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3520 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3521 if(.not.allocated(costtab)) allocate(costtab(nres))
3522 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3523 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3524 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3525 ! real(kind=8),dimension(:),allocatable :: vbld !(2*maxres) in io_conf: molread or readpdb
3526 allocate(omicron(2,nres+2)) !(2,maxres)
3527 allocate(tauangle(3,nres+2)) !(3,maxres)
3528 !elwrite(iout,*) "jestem w alloc geo 3"
3529 if(.not.allocated(xxtab)) allocate(xxtab(nres))
3530 if(.not.allocated(yytab)) allocate(yytab(nres))
3531 if(.not.allocated(zztab)) allocate(zztab(nres)) !(maxres)
3532 if(.not.allocated(xxref)) allocate(xxref(nres))
3533 if(.not.allocated(yyref)) allocate(yyref(nres))
3534 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3535 allocate(ialph(nres,2)) !(maxres,2)
3536 allocate(ivar(4*nres2)) !(4*maxres2)
3539 end subroutine alloc_geo_arrays
3540 !-----------------------------------------------------------------------------
3541 !-----------------------------------------------------------------------------