2 !-----------------------------------------------------------------------------
11 !-----------------------------------------------------------------------------
14 !-----------------------------------------------------------------------------
18 real(kind=8),dimension(:,:,:),allocatable :: t,r !(3,3,maxres)
19 !-----------------------------------------------------------------------------
22 !-----------------------------------------------------------------------------
24 ! Variables (set in init routine) never modified by local_move
26 integer :: init_called
27 logical :: locmove_output
28 real(kind=8) :: min_theta, max_theta
29 real(kind=8) :: dmin2,dmax2
30 real(kind=8) :: flag,small,small2
31 ! Workspace for local_move
33 integer :: a_n,b_n,res_n
34 real(kind=8),dimension(0:7) :: a_ang
35 real(kind=8),dimension(0:3) :: b_ang
36 real(kind=8),dimension(0:11) :: res_ang
37 logical,dimension(0:2,0:7) :: a_tab
38 logical,dimension(0:2,0:3) :: b_tab
39 logical,dimension(0:2,0:2,0:11) :: res_tab
40 !-----------------------------------------------------------------------------
41 ! integer,dimension(:),allocatable :: itype_pdb !(maxres) initialize in molread
42 !-----------------------------------------------------------------------------
45 !-----------------------------------------------------------------------------
47 !-----------------------------------------------------------------------------
49 !-----------------------------------------------------------------------------
50 real(kind=8) function ARCOS(X)
51 ! implicit real*8 (a-h,o-z)
52 ! include 'COMMON.GEO'
55 IF (DABS(X).LT.1.0D0) GOTO 1
56 ARCOS=PIPOL*(1.0d0-DSIGN(1.0D0,X))
61 !-----------------------------------------------------------------------------
63 !-----------------------------------------------------------------------------
66 ! Build the virtual polypeptide chain. Side-chain centroids are moveable.
69 ! implicit real*8 (a-h,o-z)
70 ! include 'DIMENSIONS'
71 ! include 'COMMON.CHAIN'
72 ! include 'COMMON.LOCAL'
73 ! include 'COMMON.GEO'
74 ! include 'COMMON.VAR'
75 ! include 'COMMON.IOUNITS'
76 ! include 'COMMON.NAMES'
77 ! include 'COMMON.INTERACT'
81 real(kind=8) :: be,be1,alfai
84 ! Set lprn=.true. for debugging
86 print *,"I ENTER CHAINBUILD"
88 ! Define the origin and orientation of the coordinate system and locate the
89 ! first three CA's and SC(2).
91 !elwrite(iout,*)"in chainbuild"
93 !elwrite(iout,*)"after orig_frame"
95 ! Build the alpha-carbon chain.
98 call locate_next_res(i)
100 !elwrite(iout,*)"after locate_next_res"
102 ! First and last SC must coincide with the corresponding CA.
106 dc_norm(j,nres+1)=0.0D0
107 dc(j,nres+nres)=0.0D0
108 dc_norm(j,nres+nres)=0.0D0
110 c(j,nres+nres)=c(j,nres)
113 ! Temporary diagnosis
118 write (iout,'(/a)') 'Recalculated internal coordinates'
121 c(j,nres2+2)=0.5D0*(c(j,i-1)+c(j,i+1)) !maxres2=2*maxres
124 if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i)
125 be1=rad2deg*beta(nres+i,i,nres2+2,i+1)
127 if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
128 write (iout,1212) restyp(itype(i,1),1),i,dist(i-1,i),&
129 alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,nres2+2),be1
131 1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
136 end subroutine chainbuild
137 !-----------------------------------------------------------------------------
138 subroutine orig_frame
140 ! Define the origin and orientation of the coordinate system and locate
141 ! the first three atoms.
143 ! implicit real*8 (a-h,o-z)
144 ! include 'DIMENSIONS'
145 ! include 'COMMON.CHAIN'
146 ! include 'COMMON.LOCAL'
147 ! include 'COMMON.GEO'
148 ! include 'COMMON.VAR'
151 real(kind=8) :: cost,sint
153 !el allocate(t(3,3,nres)) !(3,3,maxres)
154 !el allocate(r(3,3,nres)) !(3,3,maxres)
155 !el allocate(rt(3,3,nres)) !(3,3,maxres)
156 !el allocate(dc_norm(3,0:2*nres)) !(3,0:maxres2)
157 !el allocate(prod(3,3,nres)) !(3,3,maxres)
210 dc_norm(j,2)=prod(j,1,2)
211 dc(j,2)=vbld(3)*prod(j,1,2)
212 c(j,3)=c(j,2)+dc(j,2)
214 call locate_side_chain(2)
216 end subroutine orig_frame
217 !-----------------------------------------------------------------------------
218 subroutine locate_next_res(i)
220 ! Locate CA(i) and SC(i-1)
222 ! implicit real*8 (a-h,o-z)
223 ! include 'DIMENSIONS'
224 ! include 'COMMON.CHAIN'
225 ! include 'COMMON.LOCAL'
226 ! include 'COMMON.GEO'
227 ! include 'COMMON.VAR'
228 ! include 'COMMON.IOUNITS'
229 ! include 'COMMON.NAMES'
230 ! include 'COMMON.INTERACT'
232 ! Define the rotation matrices corresponding to CA(i)
236 real(kind=8) :: theti,phii
237 real(kind=8) :: cost,sint,cosphi,sinphi
242 call proc_proc(theti,icrc)
243 if(icrc.eq.1)theti=100.0
246 call proc_proc(phii,icrc)
247 if(icrc.eq.1)phii=180.0
250 if (theti.ne.theti) theti=100.0
252 if (phii.ne.phii) phii=180.0
262 ! Define the matrices of the rotation about the virtual-bond valence angles
263 ! theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
264 ! program), R(i,j,k), and, the cumulative matrices of rotation RT
286 rt(2,1,i-2)=sint*cosphi
287 rt(2,2,i-2)=-cost*cosphi
289 rt(3,1,i-2)=-sint*sinphi
290 rt(3,2,i-2)=cost*sinphi
292 call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
294 dc_norm(j,i-1)=prod(j,1,i-1)
295 dc(j,i-1)=vbld(i)*prod(j,1,i-1)
296 c(j,i)=c(j,i-1)+dc(j,i-1)
298 !d print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
300 ! Now calculate the coordinates of SC(i-1)
302 call locate_side_chain(i-1)
304 end subroutine locate_next_res
305 !-----------------------------------------------------------------------------
306 subroutine locate_side_chain(i)
308 ! Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.CHAIN'
313 ! include 'COMMON.LOCAL'
314 ! include 'COMMON.GEO'
315 ! include 'COMMON.VAR'
316 ! include 'COMMON.IOUNITS'
317 ! include 'COMMON.NAMES'
318 ! include 'COMMON.INTERACT'
320 real(kind=8),dimension(3) :: xx
321 real(kind=8) :: alphi,omegi,theta2
322 real(kind=8) :: dsci,dsci_inv,sinalphi,cosalphi,cosomegi,sinomegi
323 real(kind=8) :: xp,yp,zp,cost2,sint2,rj
324 ! dsci=dsc(itype(i,1))
325 ! dsci_inv=dsc_inv(itype(i,1))
327 dsci_inv=vbld_inv(i+nres)
334 call proc_proc(alphi,icrc)
335 if(icrc.eq.1)alphi=100.0
337 call proc_proc(omegi,icrc)
338 if(icrc.eq.1)omegi=-100.0
340 if (alphi.ne.alphi) alphi=100.0
341 if (omegi.ne.omegi) omegi=-100.0
352 yp= dsci*sinalphi*cosomegi
353 zp=-dsci*sinalphi*sinomegi
354 ! Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
355 ! X-axis aligned with the vector DC(*,i)
356 theta2=pi-0.5D0*theta(i+1)
359 xx(1)= xp*cost2+yp*sint2
360 xx(2)=-xp*sint2+yp*cost2
362 !d print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i,1)),i,
363 !d & xp,yp,zp,(xx(k),k=1,3)
367 ! Bring the SC vectors to the common coordinate system.
369 xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
370 xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
377 rj=rj+prod(j,k,i-1)*xx(k)
380 dc_norm(j,nres+i)=rj*dsci_inv
381 c(j,nres+i)=c(j,i)+rj
384 end subroutine locate_side_chain
385 !-----------------------------------------------------------------------------
387 !-----------------------------------------------------------------------------
388 subroutine int_from_cart1(lprn)
389 ! implicit real*8 (a-h,o-z)
390 ! include 'DIMENSIONS'
395 ! include 'COMMON.IOUNITS'
396 ! include 'COMMON.VAR'
397 ! include 'COMMON.CHAIN'
398 ! include 'COMMON.GEO'
399 ! include 'COMMON.INTERACT'
400 ! include 'COMMON.LOCAL'
401 ! include 'COMMON.NAMES'
402 ! include 'COMMON.SETUP'
403 ! include 'COMMON.TIME1'
407 real(kind=8) :: dnorm1,dnorm2,be
410 if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
417 !write(iout,*)"geometry warring, vbld=",(vbld(i),i=1,nres+1)
419 vbld_inv(nres+1)=0.0d0
420 vbld_inv(2*nres)=0.0d0
423 #if defined(PARINT) && defined(MPI)
424 do i=iint_start,iint_end
431 c(j,nres2+2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1 &
432 +(c(j,i+1)-c(j,i))/dnorm2)
436 if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
437 if ((itype(i,1).ne.10).and.(itype(i-1,1).ne.10)) then
438 tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
440 if (itype(i-1,1).ne.10) then
441 tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
442 omicron(1,i)=alpha(i-2,i-1,i-1+nres)
443 omicron(2,i)=alpha(i-1+nres,i-1,i)
445 if (itype(i,1).ne.10) then
446 tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
449 omeg(i)=beta(nres+i,i,nres2+2,i+1)
450 alph(i)=alpha(nres+i,i,nres2+2)
451 theta(i+1)=alpha(i-1,i,i+1)
453 ! print *,i,vbld(i),"vbld(i)"
454 vbld_inv(i)=1.0d0/vbld(i)
455 vbld(nres+i)=dist(nres+i,i)
456 if (itype(i,1).ne.10) then
457 vbld_inv(nres+i)=1.0d0/vbld(nres+i)
459 vbld_inv(nres+i)=0.0d0
462 #if defined(PARINT) && defined(MPI)
463 if (nfgtasks1.gt.1) then
464 !d write(iout,*) "iint_start",iint_start," iint_count",
465 !d & (iint_count(i),i=0,nfgtasks-1)," iint_displ",
466 !d & (iint_displ(i),i=0,nfgtasks-1)
467 !d write (iout,*) "Gather vbld backbone"
470 call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),&
471 MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),&
472 MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
473 !d write (iout,*) "Gather vbld_inv"
475 call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),&
476 MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),&
477 MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
478 !d write (iout,*) "Gather vbld side chain"
480 call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),&
481 MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),&
482 MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
483 !d write (iout,*) "Gather vbld_inv side chain"
485 call MPI_Allgatherv(vbld_inv(iint_start+nres),&
486 iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),&
487 iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
488 !d write (iout,*) "Gather theta"
490 call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),&
491 MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),&
492 MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
493 !d write (iout,*) "Gather phi"
495 call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),&
496 MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),&
497 MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
499 !d write (iout,*) "Gather alph"
501 call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),&
502 MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),&
503 MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
504 !d write (iout,*) "Gather omeg"
506 call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),&
507 MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),&
508 MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
510 time_gather=time_gather+MPI_Wtime()-time00
516 #if defined(WHAM_RUN) || defined(CLUSTER)
517 dc(j,i)=c(j,i+1)-c(j,i)
519 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
525 #if defined(WHAM_RUN) || defined(CLUSTER)
526 dc(j,i+nres)=c(j,i+nres)-c(j,i)
528 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
533 write (iout,1212) restyp(itype(i,1),1),i,vbld(i),&
534 rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),&
535 rad2deg*alph(i),rad2deg*omeg(i)
538 1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
540 time_intfcart=time_intfcart+MPI_Wtime()-time01
543 end subroutine int_from_cart1
544 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
545 !-----------------------------------------------------------------------------
547 !-----------------------------------------------------------------------------
548 subroutine check_sc_distr
549 ! implicit real*8 (a-h,o-z)
550 ! include 'DIMENSIONS'
551 ! include 'COMMON.TIME1'
552 ! include 'COMMON.INTERACT'
553 ! include 'COMMON.NAMES'
554 ! include 'COMMON.GEO'
555 ! include 'COMMON.HEADER'
556 ! include 'COMMON.CONTROL'
558 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
559 real(kind=8) :: hrtime,mintime,sectime
560 integer,parameter :: MaxSample=10000000
561 real(kind=8),parameter :: delt=1.0D0/MaxSample
562 real(kind=8),dimension(0:72,0:90) :: prob
564 integer :: it,i,j,isample,indal,indom
565 real(kind=8) :: al,om,dV
566 dV=2.0D0*5.0D0*deg2rad*deg2rad
569 if ((it.eq.10).or.(it.eq.ntyp1)) goto 10
570 open (20,file=restyp(it,1)//'_distr.sdc',status='unknown')
571 call gen_side(it,90.0D0 * deg2rad,al,om,fail,1)
574 open (20,file=restyp(it,1)//'_distr1.sdc',status='unknown')
580 do isample=1,MaxSample
581 call gen_side(it,90.0D0 * deg2rad,al,om,fail,1)
583 indom=(rad2deg*om+180.0D0)/5
584 prob(indom,indal)=prob(indom,indal)+delt
588 write (20,'(2f10.3,1pd15.5)') 2*i+0.0D0,5*j-180.0D0,&
594 end subroutine check_sc_distr
596 !-----------------------------------------------------------------------------
598 !-----------------------------------------------------------------------------
599 subroutine geom_to_var(n,x)
601 ! Transfer the geometry parameters to the variable array.
602 ! The positions of variables are as follows:
603 ! 1. Virtual-bond torsional angles: 1 thru nres-3
604 ! 2. Virtual-bond valence angles: nres-2 thru 2*nres-5
605 ! 3. The polar angles alpha of local SC orientation: 2*nres-4 thru
607 ! 4. The torsional angles omega of SC orientation: 2*nres-4+nside+1
608 ! thru 2*nre-4+2*nside
610 ! implicit real*8 (a-h,o-z)
611 ! include 'DIMENSIONS'
612 ! include 'COMMON.VAR'
613 ! include 'COMMON.GEO'
614 ! include 'COMMON.CHAIN'
616 real(kind=8),dimension(n) :: x
617 !d print *,'nres',nres,' nphi',nphi,' ntheta',ntheta,' nvar',nvar
620 !d print *,i,i-3,phi(i)
622 if (n.eq.nphi) return
625 !d print *,i,i-2+nphi,theta(i)
627 if (n.eq.nphi+ntheta) return
629 if (ialph(i,1).gt.0) then
630 x(ialph(i,1))=alph(i)
631 x(ialph(i,1)+nside)=omeg(i)
632 !d print *,i,ialph(i,1),ialph(i,1)+nside,alph(i),omeg(i)
636 end subroutine geom_to_var
637 !-----------------------------------------------------------------------------
638 subroutine var_to_geom(n,x)
640 ! Update geometry parameters according to the variable array.
642 ! implicit real*8 (a-h,o-z)
643 ! include 'DIMENSIONS'
644 ! include 'COMMON.VAR'
645 ! include 'COMMON.CHAIN'
646 ! include 'COMMON.GEO'
647 ! include 'COMMON.IOUNITS'
649 real(kind=8),dimension(n) :: x
650 logical :: change !,reduce
657 if (n.gt.nphi+ntheta) then
660 alph(ii)=x(nphi+ntheta+i)
661 omeg(ii)=pinorm(x(nphi+ntheta+nside+i))
662 !elwrite(iout,*) "alph",ii,alph
663 !elwrite(iout,*) "omeg",ii,omeg
668 !elwrite(iout,*) "phi",i,phi
670 if (n.eq.nphi) return
673 !elwrite(iout,*) "theta",i,theta
674 if (theta(i).eq.pi) theta(i)=0.99d0*pi
678 end subroutine var_to_geom
679 !-----------------------------------------------------------------------------
680 logical function convert_side(alphi,omegi)
682 real(kind=8) :: alphi,omegi
683 !el real(kind=8) :: pinorm
684 ! include 'COMMON.GEO'
686 ! Apply periodicity restrictions.
687 if (alphi.gt.pi) then
689 omegi=pinorm(omegi+pi)
693 end function convert_side
694 !-----------------------------------------------------------------------------
695 logical function reduce(x)
697 ! Apply periodic restrictions to variables.
699 ! implicit real*8 (a-h,o-z)
700 ! include 'DIMENSIONS'
701 ! include 'COMMON.VAR'
702 ! include 'COMMON.CHAIN'
703 ! include 'COMMON.GEO'
704 logical :: zm,zmiana !,convert_side
705 real(kind=8),dimension(nvar) :: x
709 x(i-3)=pinorm(x(i-3))
711 if (nvar.gt.nphi+ntheta) then
715 x(ii)=thetnorm(x(ii))
716 x(iii)=pinorm(x(iii))
717 ! Apply periodic restrictions.
718 zm=convert_side(x(ii),x(iii))
722 if (nvar.eq.nphi) return
726 x(ii)=dmod(x(ii),dwapi)
727 ! Apply periodic restrictions.
728 if (x(ii).gt.pi) then
731 if (iii.gt.0) x(iii)=pinorm(x(iii)+pi)
732 if (i.lt.nres) x(iii+1)=pinorm(x(iii+1)+pi)
735 x(ii)=dmod(pi-x(ii),dwapi)
736 x(ii+nside)=pinorm(-x(ii+nside))
737 zm=convert_side(x(ii),x(ii+nside))
739 else if (x(ii).lt.-pi) then
744 x(ii)=dmod(pi-x(ii),dwapi)
745 x(ii+nside)=pinorm(-pi-x(ii+nside))
746 zm=convert_side(x(ii),x(ii+nside))
748 else if (x(ii).lt.0.0d0) then
751 if (iii.gt.0) x(iii)=pinorm(x(iii)+pi)
752 if (i.lt.nres) x(iii+1)=pinorm(x(iii+1)+pi)
755 x(ii+nside)=pinorm(-x(ii+nside))
756 zm=convert_side(x(ii),x(ii+nside))
763 !-----------------------------------------------------------------------------
764 real(kind=8) function thetnorm(x)
765 ! This function puts x within [0,2Pi].
768 ! include 'COMMON.GEO'
770 if (xx.lt.0.0d0) xx=xx+dwapi
771 if (xx.gt.0.9999d0*pi) xx=0.9999d0*pi
774 end function thetnorm
775 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
776 !-----------------------------------------------------------------------------
777 subroutine var_to_geom_restr(n,xx)
779 ! Update geometry parameters according to the variable array.
781 ! implicit real*8 (a-h,o-z)
782 ! include 'DIMENSIONS'
783 ! include 'COMMON.VAR'
784 ! include 'COMMON.CHAIN'
785 ! include 'COMMON.GEO'
786 ! include 'COMMON.IOUNITS'
788 real(kind=8),dimension(6*nres) :: x,xx !(maxvar) (maxvar=6*maxres)
789 logical :: change !,reduce
795 alph(ii)=x(nphi+ntheta+i)
796 omeg(ii)=pinorm(x(nphi+ntheta+nside+i))
803 if (theta(i).eq.pi) theta(i)=0.99d0*pi
807 end subroutine var_to_geom_restr
808 !-----------------------------------------------------------------------------
810 !-----------------------------------------------------------------------------
811 subroutine gen_rand_conf(nstart,*)
812 ! Generate random conformation or chain cut and regrowth.
814 use random, only: iran_num,ran_number
815 ! implicit real*8 (a-h,o-z)
816 ! include 'DIMENSIONS'
817 ! include 'COMMON.CHAIN'
818 ! include 'COMMON.LOCAL'
819 ! include 'COMMON.VAR'
820 ! include 'COMMON.INTERACT'
821 ! include 'COMMON.IOUNITS'
822 ! include 'COMMON.MCM'
823 ! include 'COMMON.GEO'
824 ! include 'COMMON.CONTROL'
825 logical :: back,fail !overlap,
827 integer :: i,nstart,maxsi,nsi,maxnit,nit,niter
828 integer :: it1,it2,it,j
829 !d print *,' CG Processor',me,' maxgen=',maxgen
831 write (iout,*) 'Gen_Rand_conf: nstart=',nstart,nres
832 if (nstart.lt.5) then
834 phi(4)=gen_phi(4,iabs(itype(2,1)),iabs(itype(3,1)))
835 ! write(iout,*)'phi(4)=',rad2deg*phi(4)
836 if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2,1)),pi,phi(4),molnum(2))
837 ! write(iout,*)'theta(3)=',rad2deg*theta(3)
838 if ((it1.ne.10).and.(it1.ne.ntyp1)) then
841 do while (fail.and.nsi.le.maxsi)
842 call gen_side(it1,theta(3),alph(2),omeg(2),fail,molnum(2))
843 write (iout,*) 'nsi=',nsi,maxsi
846 if (nsi.gt.maxsi) return 1
848 write(iout,*) "before origin_frame"
850 write(iout,*) "after origin_frame"
863 do while (i.le.nres .and. niter.lt.maxgen)
864 write(iout,*) 'i=',i,'back=',back
865 if (i.lt.nstart) then
867 write (iout,'(/80(1h*)/2a/80(1h*))') &
868 'Generation procedure went down to ',&
869 'chain beginning. Cannot continue...'
870 write (*,'(/80(1h*)/2a/80(1h*))') &
871 'Generation procedure went down to ',&
872 'chain beginning. Cannot continue...'
876 it1=iabs(itype(i-1,molnum(i-1)))
877 it2=iabs(itype(i-2,molnum(i-2)))
878 it=iabs(itype(i,molnum(i)))
879 if ((it.eq.ntyp1).and.(it1.eq.ntyp1)) &
880 vbld(i)=ran_number(30.0D0,40.0D0)
881 ! print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,&
882 ! ' nit=',nit,' niter=',niter,' maxgen=',maxgen
883 phi(i+1)=gen_phi(i+1,it1,it)
885 phi(i)=gen_phi(i+1,it2,it1)
886 ! print *,'phi(',i,')=',phi(i)
887 theta(i-1)=gen_theta(it2,phi(i-1),phi(i),molnum(i))
888 ! print *,"theta",theta(i-1),phi(i)
889 if ((it2.ne.10).and.(it2.ne.ntyp1)) then
892 do while (fail.and.nsi.le.maxsi)
893 call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail,molnum(i-2))
896 if (nsi.gt.maxsi) return 1
898 call locate_next_res(i-1)
900 theta(i)=gen_theta(it1,phi(i),phi(i+1),molnum(i))
901 ! write(iout,*) "theta(i),",theta(i)
902 if ((it1.ne.10).and.(it1.ne.ntyp1)) then
905 do while (fail.and.nsi.le.maxsi)
906 call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail,molnum(i))
907 ! write(iout,*)"alpha,omeg(i-1)",alph(i-1),omeg(i-1),i,nsi,maxsi
910 if (nsi.gt.maxsi) return 1
912 call locate_next_res(i)
913 write(iout,*) "overlap,",overlap(i-1)
914 if (overlap(i-1)) then
915 if (nit.lt.maxnit) then
925 'Cannot generate non-overlaping conformation. Increase MAXNIT.'
927 'Cannot generate non-overlaping conformation. Increase MAXNIT.'
932 ! write(iout,*) "tu dochodze"
939 if (niter.ge.maxgen) then
940 write (iout,'(a,2i5)') &
941 'Too many trials in conformation generation',niter,maxgen
942 write (*,'(a,2i5)') &
943 'Too many trials in conformation generation',niter,maxgen
948 c(j,nres+nres)=c(j,nres)
951 end subroutine gen_rand_conf
952 !-----------------------------------------------------------------------------
953 logical function overlap(i)
954 ! implicit real*8 (a-h,o-z)
955 ! include 'DIMENSIONS'
956 ! include 'COMMON.CHAIN'
957 ! include 'COMMON.INTERACT'
958 ! include 'COMMON.FFIELD'
959 integer :: i,j,iti,itj,iteli,itelj,k
960 real(kind=8) :: redfac,rcomp
965 iti=iabs(itype(i,molnum(i)))
966 if (iti.gt.ntyp) return
967 ! Check for SC-SC overlaps.
968 !d print *,'nnt=',nnt,' nct=',nct
970 ! print *, "molnum(j)",j,molnum(j)
971 if (molnum(j).eq.1) then
973 if (itj.eq.ntyp1) cycle
974 if (j.lt.i-1 .or. ipot.ne.4) then
975 rcomp=sigmaii(iti,itj)
980 if (dist(nres+i,nres+j).lt.redfac*rcomp) then
983 ! print *,'overlap, SC-SC: i=',i,' j=',j,
984 ! & ' dist=',dist(nres+i,nres+j),' rcomp=',
988 else if (molnum(j).eq.2) then
990 if (dist(nres+i,nres+j).lt.redfac*sigma_nucl(iti,itj)) then
993 ! print *,'overlap, SC-SC: i=',i,' j=',j,
994 ! & ' dist=',dist(nres+i,nres+j),' rcomp=',
1001 ! Check for overlaps between the added peptide group and the preceding
1005 ! c(j,nres2+1)=0.5D0*(c(j,i)+c(j,i+1))
1006 c(j,nres2+3)=0.5D0*(c(j,i)+c(j,i+1))
1009 if (molnum(j).ne.1) cycle
1010 itj=iabs(itype(j,1))
1011 !d print *,'overlap, p-Sc: i=',i,' j=',j,
1012 !d & ' dist=',dist(nres+j,maxres2+1)
1013 if (dist(nres+j,nres2+3).lt.4.0D0*redfac) then
1018 ! Check for overlaps between the added side chain and the preceding peptide
1021 if (molnum(j).ne.1) cycle
1023 c(k,nres2+3)=0.5D0*(c(k,j)+c(k,j+1))
1025 !d print *,'overlap, SC-p: i=',i,' j=',j,
1026 !d & ' dist=',dist(nres+i,maxres2+1)
1027 if (dist(nres+i,nres2+3).lt.4.0D0*redfac) then
1032 ! Check for p-p overlaps
1034 c(j,nres2+3)=0.5D0*(c(j,i)+c(j,i+1))
1037 ! if (molnum(j).eq.1) then
1040 c(k,nres2+4)=0.5D0*(c(k,j)+c(k,j+1))
1042 !d print *,'overlap, p-p: i=',i,' j=',j,
1043 !d & ' dist=',dist(maxres2+1,maxres2+2)
1044 if (molnum(j).eq.1) then
1045 if(iteli.ne.0.and.itelj.ne.0)then
1046 if (dist(nres2+3,nres2+4).lt.rpp(iteli,itelj)*redfac) then
1051 else if (molnum(j).eq.2) then
1052 if (dist(nres2+3,nres2+4).lt.3.0) then
1059 end function overlap
1060 !-----------------------------------------------------------------------------
1061 real(kind=8) function gen_phi(i,it1,it2)
1062 use random, only:ran_number
1063 ! implicit real*8 (a-h,o-z)
1064 ! include 'DIMENSIONS'
1065 ! include 'COMMON.GEO'
1066 ! include 'COMMON.BOUNDS'
1067 integer :: i,it1,it2
1068 ! gen_phi=ran_number(-pi,pi)
1069 ! 8/13/98 Generate phi using pre-defined boundaries
1070 gen_phi=ran_number(phibound(1,i),phibound(2,i))
1072 end function gen_phi
1073 !-----------------------------------------------------------------------------
1074 real(kind=8) function gen_theta(it,gama,gama1,mnum)
1075 use random,only:binorm,ran_number
1076 ! implicit real*8 (a-h,o-z)
1077 ! include 'DIMENSIONS'
1078 ! include 'COMMON.LOCAL'
1079 ! include 'COMMON.GEO'
1080 real(kind=8),dimension(2) :: y,z
1081 real(kind=8) :: theta_max,theta_min,sig,ak
1083 integer :: j,it,k,mnum
1084 real(kind=8) :: gama,gama1,thet_pred_mean,theta_temp
1085 ! print *,'gen_theta: it=',it
1088 if (dabs(gama).gt.dwapi) then
1095 if (dabs(gama1).gt.dwapi) then
1102 if (it.eq.ntyp1) then
1103 gen_theta=ran_number(theta_max/2.0,theta_max)
1104 else if (mnum.eq.1) then
1106 thet_pred_mean=a0thet(it)
1107 ! write(iout,*),it,thet_pred_mean,"gen_thet"
1109 thet_pred_mean=thet_pred_mean+athet(k,it,1,1)*y(k) &
1110 +bthet(k,it,1,1)*z(k)
1114 sig=sig*thet_pred_mean+polthet(j,it)
1116 sig=0.5D0/(sig*sig+sigc0(it))
1117 ak=dexp(gthet(1,it)- &
1118 0.5D0*((gthet(2,it)-thet_pred_mean)/gthet(3,it))**2)
1119 ! print '(i5,5(1pe14.4))',it,(gthet(j,it),j=1,3)
1120 ! print '(5(1pe14.4))',thet_pred_mean,theta0(it),sig,sig0(it),ak
1121 theta_temp=binorm(thet_pred_mean,theta0(it),sig,sig0(it),ak)
1122 if (theta_temp.lt.theta_min) theta_temp=theta_min
1123 if (theta_temp.gt.theta_max) theta_temp=theta_max
1124 gen_theta=theta_temp
1125 ! print '(a)','Exiting GENTHETA.'
1126 else if (mnum.eq.2) then
1127 gen_theta=2.0d0 + ran_number(0.0d0,0.34d0)
1129 gen_theta=ran_number(theta_max/2.0,theta_max)
1132 end function gen_theta
1133 !-----------------------------------------------------------------------------
1134 subroutine gen_side(it,the,al,om,fail,mnum)
1135 use random, only:ran_number,mult_norm1
1136 ! implicit real*8 (a-h,o-z)
1137 ! include 'DIMENSIONS'
1138 ! include 'COMMON.GEO'
1139 ! include 'COMMON.LOCAL'
1140 ! include 'COMMON.SETUP'
1141 ! include 'COMMON.IOUNITS'
1142 real(kind=8) :: MaxBoxLen=10.0D0
1143 real(kind=8),dimension(3,3) :: Ap_inv,a,vec
1144 real(kind=8),dimension(:,:),allocatable :: z !(3,maxlob)
1145 real(kind=8),dimension(:),allocatable :: W1,detAp !(maxlob)
1146 real(kind=8),dimension(:),allocatable :: sumW !(0:maxlob)
1147 real(kind=8),dimension(2) :: y,cm,eig
1148 real(kind=8),dimension(2,2) :: box
1149 real(kind=8),dimension(100) :: work
1150 real(kind=8) :: eig_limit=1.0D-8
1151 real(kind=8) :: Big=10.0D0
1152 logical :: lprint,fail,lcheck
1154 integer :: it,i,j,k,l,nlobit,ial,iom,iii,ilob,mnum
1155 real(kind=8) :: the,al,om,detApi,wart,y2,wykl,radmax
1156 real(kind=8) :: tant,zz1,W1i,radius,zk,fac,dV,sum,sum1
1157 real(kind=8) :: which_lobe
1162 if (the.eq.0.0D0 .or. the.eq.pi) then
1164 write (*,'(a,i4,a,i3,a,1pe14.5)') &
1165 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the
1167 !d write (iout,'(a,i3,a,1pe14.5)')
1168 !d & 'Error in GenSide: it=',it,' theta=',the
1173 if (nlobit.eq.0) then
1174 al=ran_number(0.05d0,pi/2)
1175 om=ran_number(-pi,pi)
1178 tant=dtan(the-pipol)
1180 allocate(z(3,nlobit))
1181 allocate(W1(nlobit))
1182 allocate(detAp(nlobit))
1183 allocate(sumW(0:nlobit))
1186 print '(a,i4,a)','CG Processor:',me,' Enter Gen_Side.'
1187 write (iout,'(a,i4,a)') 'Processor:',me,' Enter Gen_Side.'
1189 print *,'it=',it,' nlobit=',nlobit,' the=',the,' tant=',tant
1190 write (iout,*) 'it=',it,' nlobit=',nlobit,' the=',the,&
1194 zz1=tant-censc(1,i,it)
1197 a(k,l)=gaussc(k,l,i,it)
1200 detApi=a(2,2)*a(3,3)-a(2,3)**2
1201 Ap_inv(2,2)=a(3,3)/detApi
1202 Ap_inv(2,3)=-a(2,3)/detApi
1203 Ap_inv(3,2)=Ap_inv(2,3)
1204 Ap_inv(3,3)=a(2,2)/detApi
1206 write (*,'(/a,i2/)') 'Cluster #',i
1207 write (*,'(3(1pe14.5),5x,1pe14.5)') &
1208 ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
1209 write (iout,'(/a,i2/)') 'Cluster #',i
1210 write (iout,'(3(1pe14.5),5x,1pe14.5)') &
1211 ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
1216 W1i=W1i+a(k,1)*a(l,1)*Ap_inv(k,l)
1220 W1(i)=dexp(bsc(i,it)-0.5D0*W1i*zz1*zz1)
1221 ! if (lprint) write(*,'(a,3(1pe15.5)/)')
1222 ! & 'detAp, W1, anormi',detApi,W1i,anormi
1226 zk=zk+zz1*Ap_inv(k,l)*a(l,1)
1230 detAp(i)=dsqrt(detApi)
1234 print *,'W1:',(w1(i),i=1,nlobit)
1235 print *,'detAp:',(detAp(i),i=1,nlobit)
1238 print '(i2,3f10.5)',i,(rad2deg*z(j,i),j=2,3)
1240 write (iout,*) 'W1:',(w1(i),i=1,nlobit)
1241 write (iout,*) 'detAp:',(detAp(i),i=1,nlobit)
1244 write (iout,'(i2,3f10.5)') i,(rad2deg*z(j,i),j=2,3)
1248 ! Writing the distribution just to check the procedure
1250 dV=deg2rad**2*10.0D0
1254 fac=fac+W1(i)/detAp(i)
1256 fac=1.0D0/(2.0D0*fac*pi)
1257 !d print *,it,'fac=',fac
1266 a(j-1,k-1)=gaussc(j,k,i,it)
1278 wykl=wykl+a(j,k)*(y(j)-z(j+1,i))*(y(k)-z(k+1,i))
1281 wart=wart+W1(i)*dexp(-0.5D0*wykl)
1288 ! print *,'y',y(1),y(2),' fac=',fac
1290 write (20,'(2f10.3,1pd15.5)') y(1)*rad2deg,y(2)*rad2deg,wart
1295 ! print *,'it=',it,' sum=',sum*dV,' sum1=',sum1*dV
1299 ! Calculate the CM of the system
1302 W1(i)=W1(i)/detAp(i)
1306 sumW(i)=sumW(i-1)+W1(i)
1311 cm(1)=cm(1)+z(2,j)*W1(j)
1312 cm(2)=cm(2)+W1(j)*(z(3,1)+pinorm(z(3,j)-z(3,1)))
1314 cm(1)=cm(1)/sumW(nlobit)
1315 cm(2)=cm(2)/sumW(nlobit)
1316 if (cm(1).gt.Big .or. cm(1).lt.-Big .or. &
1317 cm(2).gt.Big .or. cm(2).lt.-Big) then
1318 !d write (iout,'(a)')
1319 !d & 'Unexpected error in GenSide - CM coordinates too large.'
1320 !d write (iout,'(i5,2(1pe14.5))') it,cm(1),cm(2)
1322 !d & 'Unexpected error in GenSide - CM coordinates too large.'
1323 !d write (*,'(i5,2(1pe14.5))') it,cm(1),cm(2)
1327 !d print *,'CM:',cm(1),cm(2)
1329 ! Find the largest search distance from CM
1335 a(j-1,k-1)=gaussc(j,k,i,it)
1339 call f02faf('N','U',2,a,3,eig,work,100,ifail)
1341 call djacob(2,3,10000,1.0d-10,a,vec,eig)
1345 print *,'*************** CG Processor',me
1346 print *,'CM:',cm(1),cm(2)
1347 write (iout,*) '*************** CG Processor',me
1348 write (iout,*) 'CM:',cm(1),cm(2)
1349 print '(A,8f10.5)','Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
1350 write (iout,'(A,8f10.5)') &
1351 'Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
1354 if (eig(1).lt.eig_limit) then
1356 'From Mult_Norm: Eigenvalues of A are too small.'
1358 'From Mult_Norm: Eigenvalues of A are too small.'
1365 radius=radius+pinorm(z(j+1,i)-cm(j))**2
1367 radius=dsqrt(radius)+3.0D0/dsqrt(eig(1))
1368 if (radius.gt.radmax) radmax=radius
1370 if (radmax.gt.pi) radmax=pi
1372 ! Determine the boundaries of the search rectangle.
1375 print '(a,4(1pe14.4))','W1: ',(W1(i),i=1,nlob(it) )
1376 print '(a,4(1pe14.4))','radmax: ',radmax
1378 box(1,1)=dmax1(cm(1)-radmax,0.0D0)
1379 box(2,1)=dmin1(cm(1)+radmax,pi)
1380 box(1,2)=cm(2)-radmax
1381 box(2,2)=cm(2)+radmax
1384 print *,'CG Processor',me,' Array BOX:'
1386 print *,'Array BOX:'
1388 print '(4(1pe14.4))',((box(k,j),k=1,2),j=1,2)
1389 print '(a,4(1pe14.4))','sumW: ',(sumW(i),i=0,nlob(it) )
1391 write (iout,*)'CG Processor',me,' Array BOX:'
1393 write (iout,*)'Array BOX:'
1395 write(iout,'(4(1pe14.4))') ((box(k,j),k=1,2),j=1,2)
1396 write(iout,'(a,4(1pe14.4))')'sumW: ',(sumW(i),i=0,nlob(it) )
1398 ! if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
1400 ! write (iout,'(a,i4,a,3e15.5)') 'CG Processor:',me,': bad sampling box.',box(1,2),box(2,2),radmax
1401 ! write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
1403 ! write (iout,'(a)') 'Bad sampling box.'
1408 which_lobe=ran_number(0.0D0,sumW(nlobit))
1409 ! print '(a,1pe14.4)','which_lobe=',which_lobe
1411 if (sumW(i-1).le.which_lobe .and. sumW(i).ge.which_lobe) goto 1
1414 ! print *,'ilob=',ilob,' nlob=',nlob(it)
1418 a(i-1,j-1)=gaussc(i,j,ilob,it)
1421 !d print '(a,i4,a)','CG Processor',me,' Calling MultNorm1.'
1422 call mult_norm1(3,2,a,cm,box,y,fail)
1426 else if (mnum.eq.2) then
1427 al=0.7+ran_number(0.0d0,0.2d0)
1428 om=ran_number(0.0d0,3.14d0)
1431 !d print *,'al=',al,' om=',om
1434 end subroutine gen_side
1435 !-----------------------------------------------------------------------------
1436 subroutine overlap_sc(scfail)
1438 ! Internal and cartesian coordinates must be consistent as input,
1439 ! and will be up-to-date on return.
1440 ! At the end of this procedure, scfail is true if there are
1441 ! overlapping residues left, or false otherwise (success)
1443 ! implicit real*8 (a-h,o-z)
1444 ! include 'DIMENSIONS'
1445 ! include 'COMMON.CHAIN'
1446 ! include 'COMMON.INTERACT'
1447 ! include 'COMMON.FFIELD'
1448 ! include 'COMMON.VAR'
1449 ! include 'COMMON.SBRIDGE'
1450 ! include 'COMMON.IOUNITS'
1451 logical :: had_overlaps,fail,scfail
1452 integer,dimension(nres) :: ioverlap !(maxres)
1453 integer :: ioverlap_last,k,maxsi,i,iti,nsi
1456 had_overlaps=.false.
1457 call overlap_sc_list(ioverlap,ioverlap_last)
1458 if (ioverlap_last.gt.0) then
1459 write (iout,*) '#OVERLAPing residues ',ioverlap_last
1460 write (iout,'(20i4)') (ioverlap(k),k=1,ioverlap_last)
1466 if (ioverlap_last.eq.0) exit
1468 do ires=1,ioverlap_last
1470 iti=iabs(itype(i,1))
1471 if ((iti.ne.10).and.(molnum(i).ne.5).and.(iti.ne.ntyp1)) then
1474 do while (fail.and.nsi.le.maxsi)
1475 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail,molnum(i))
1483 call overlap_sc_list(ioverlap,ioverlap_last)
1484 ! write (iout,*) 'Overlaping residues ',ioverlap_last,
1485 ! & (ioverlap(j),j=1,ioverlap_last)
1488 if (k.le.1000.and.ioverlap_last.eq.0) then
1490 if (had_overlaps) then
1491 write (iout,*) '#OVERLAPing all corrected after ',k,&
1492 ' random generation'
1496 write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
1497 write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last)
1503 write (iout,'(a30,i5,a12,i4)') &
1504 '#OVERLAP FAIL in gen_side after',maxsi,&
1508 end subroutine overlap_sc
1509 !-----------------------------------------------------------------------------
1510 subroutine overlap_sc_list(ioverlap,ioverlap_last)
1512 ! implicit real*8 (a-h,o-z)
1513 ! include 'DIMENSIONS'
1514 ! include 'COMMON.GEO'
1515 ! include 'COMMON.LOCAL'
1516 ! include 'COMMON.IOUNITS'
1517 ! include 'COMMON.CHAIN'
1518 ! include 'COMMON.INTERACT'
1519 ! include 'COMMON.FFIELD'
1520 ! include 'COMMON.VAR'
1521 ! include 'COMMON.CALC'
1523 integer,dimension(nres) :: ioverlap !(maxres)
1524 integer :: ioverlap_last
1527 real(kind=8) :: redfac,sig !rrij,sigsq,
1528 integer :: itypi,itypj,itypi1
1529 real(kind=8) :: xi,yi,zi,sig0ij,rcomp,rrij,rij_shift
1533 ! Check for SC-SC overlaps and mark residues
1534 ! print *,'>>overlap_sc nnt=',nnt,' nct=',nct
1536 do i=iatsc_s,iatsc_e
1537 if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) cycle
1538 if (molnum(i).eq.5) print *,"WTF",i,iatsc_s,iatsc_e
1539 if (molnum(i).eq.5) cycle
1540 itypi=iabs(itype(i,molnum(i)))
1541 itypi1=iabs(itype(i+1,1))
1545 dxi=dc_norm(1,nres+i)
1546 dyi=dc_norm(2,nres+i)
1547 dzi=dc_norm(3,nres+i)
1548 dsci_inv=dsc_inv(itypi)
1550 do iint=1,nint_gr(i)
1551 do j=istart(i,iint),iend(i,iint)
1552 if (itype(j,molnum(j)).eq.ntyp1_molec(molnum(j))) cycle
1554 itypj=iabs(itype(j,molnum(j)))
1555 dscj_inv=dsc_inv(itypj)
1556 sig0ij=sigma(itypi,itypj)
1557 chi1=chi(itypi,itypj)
1558 chi2=chi(itypj,itypi)
1565 alf12=0.5D0*(alf1+alf2)
1567 rcomp=sigmaii(itypi,itypj)
1569 rcomp=sigma(itypi,itypj)
1571 ! print '(2(a3,2i3),a3,2f10.5)',
1572 ! & ' i=',i,iti,' j=',j,itj,' d=',dist(nres+i,nres+j)
1577 dxj=dc_norm(1,nres+j)
1578 dyj=dc_norm(2,nres+j)
1579 dzj=dc_norm(3,nres+j)
1580 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1584 sig=sig0ij*dsqrt(sigsq)
1585 rij_shift=1.0D0/rij-sig+sig0ij
1587 !t if ( 1.0/rij .lt. redfac*rcomp .or.
1588 !t & rij_shift.le.0.0D0 ) then
1589 if ( rij_shift.le.0.0D0 ) then
1590 !d write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)')
1591 !d & 'overlap SC-SC: i=',i,' j=',j,
1592 !d & ' dist=',dist(nres+i,nres+j),' rcomp=',
1593 !d & rcomp,1.0/rij,rij_shift
1594 ioverlap_last=ioverlap_last+1
1595 ioverlap(ioverlap_last)=i
1596 do k=1,ioverlap_last-1
1597 if (ioverlap(k).eq.i) ioverlap_last=ioverlap_last-1
1599 ioverlap_last=ioverlap_last+1
1600 ioverlap(ioverlap_last)=j
1601 do k=1,ioverlap_last-1
1602 if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1
1609 end subroutine overlap_sc_list
1611 !-----------------------------------------------------------------------------
1612 ! energy_p_new_barrier.F
1613 !-----------------------------------------------------------------------------
1614 subroutine sc_angular
1615 ! Calculate eps1,eps2,eps3,sigma, and parts of their derivatives in om1,om2,
1616 ! om12. Called by ebp, egb, and egbv.
1619 ! include 'COMMON.CALC'
1620 ! include 'COMMON.IOUNITS'
1624 om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
1625 om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
1626 om12=dxi*dxj+dyi*dyj+dzi*dzj
1628 ! Calculate eps1(om12) and its derivative in om12
1629 faceps1=1.0D0-om12*chiom12
1630 faceps1_inv=1.0D0/faceps1
1631 eps1=dsqrt(faceps1_inv)
1632 ! Following variable is eps1*deps1/dom12
1633 eps1_om12=faceps1_inv*chiom12
1638 ! write (iout,*) "om12",om12," eps1",eps1
1639 ! Calculate sigma(om1,om2,om12) and the derivatives of sigma**2 in om1,om2,
1644 facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
1645 ! print *,"TUT?",om1*chiom1,facsig,om1,om2,om12
1646 sigsq=1.0D0-facsig*faceps1_inv
1647 sigsq_om1=(chiom1-chiom12*om2)*faceps1_inv
1648 sigsq_om2=(chiom2-chiom12*om1)*faceps1_inv
1649 sigsq_om12=-chi12*(om1om2*faceps1-om12*facsig)*faceps1_inv**2
1655 ! write (iout,*) "chiom1",chiom1," chiom2",chiom2," chiom12",chiom12
1656 ! write (iout,*) "faceps1",faceps1," faceps1_inv",faceps1_inv,
1658 ! Calculate eps2 and its derivatives in om1, om2, and om12.
1661 chipom12=chip12*om12
1662 facp=1.0D0-om12*chipom12
1664 facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12
1665 ! write (iout,*) "chipom1",chipom1," chipom2",chipom2,
1666 ! & " chipom12",chipom12," facp",facp," facp_inv",facp_inv
1667 ! Following variable is the square root of eps2
1668 eps2rt=1.0D0-facp1*facp_inv
1669 ! Following three variables are the derivatives of the square root of eps
1670 ! in om1, om2, and om12.
1671 eps2rt_om1=-4.0D0*(chipom1-chipom12*om2)*facp_inv
1672 eps2rt_om2=-4.0D0*(chipom2-chipom12*om1)*facp_inv
1673 eps2rt_om12=4.0D0*chip12*(om1om2*facp-om12*facp1)*facp_inv**2
1674 ! Evaluate the "asymmetric" factor in the VDW constant, eps3
1675 eps3rt=1.0D0-alf1*om1+alf2*om2-alf12*om12
1676 ! write (iout,*) "eps2rt",eps2rt," eps3rt",eps3rt
1677 ! write (iout,*) "eps2rt_om1",eps2rt_om1," eps2rt_om2",eps2rt_om2,
1678 ! & " eps2rt_om12",eps2rt_om12
1679 ! Calculate whole angle-dependent part of epsilon and contributions
1680 ! to its derivatives
1682 end subroutine sc_angular
1683 !-----------------------------------------------------------------------------
1685 subroutine sc_angular_nucl
1686 ! Calculate eps1,eps2,eps3,sigma, and parts of their derivatives in om1,om2,
1687 ! om12. Called by ebp, egb, and egbv.
1690 ! include 'COMMON.CALC'
1691 ! include 'COMMON.IOUNITS'
1697 om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
1698 om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
1699 om12=dxi*dxj+dyi*dyj+dzi*dzj
1701 ! Calculate eps1(om12) and its derivative in om12
1702 faceps1=1.0D0-om12*chiom12
1703 faceps1_inv=1.0D0/faceps1
1704 eps1=dsqrt(faceps1_inv)
1705 ! Following variable is eps1*deps1/dom12
1706 eps1_om12=faceps1_inv*chiom12
1711 ! write (iout,*) "om12",om12," eps1",eps1
1712 ! Calculate sigma(om1,om2,om12) and the derivatives of sigma**2 in om1,om2,
1717 facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
1718 sigsq=1.0D0-facsig*faceps1_inv
1719 sigsq_om1=(chiom1-chiom12*om2)*faceps1_inv
1720 sigsq_om2=(chiom2-chiom12*om1)*faceps1_inv
1721 sigsq_om12=-chi12*(om1om2*faceps1-om12*facsig)*faceps1_inv**2
1724 chipom12=chip12*om12
1725 facp=1.0D0-om12*chipom12
1727 facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12
1728 ! write (iout,*) "chipom1",chipom1," chipom2",chipom2,
1729 ! & " chipom12",chipom12," facp",facp," facp_inv",facp_inv
1730 ! Following variable is the square root of eps2
1731 eps2rt=1.0D0-facp1*facp_inv
1732 ! Following three variables are the derivatives of the square root of eps
1733 ! in om1, om2, and om12.
1734 eps2rt_om1=-4.0D0*(chipom1-chipom12*om2)*facp_inv
1735 eps2rt_om2=-4.0D0*(chipom2-chipom12*om1)*facp_inv
1736 eps2rt_om12=4.0D0*chip12*(om1om2*facp-om12*facp1)*facp_inv**2
1737 ! Evaluate the "asymmetric" factor in the VDW constant, eps3
1738 eps3rt=1.0D0-alf1*om1+alf2*om2-alf12*om12
1739 ! write (iout,*) "eps2rt",eps2rt," eps3rt",eps3rt
1740 ! write (iout,*) "eps2rt_om1",eps2rt_om1," eps2rt_om2",eps2rt_om2,
1741 ! & " eps2rt_om12",eps2rt_om12
1742 ! Calculate whole angle-dependent part of epsilon and contributions
1743 ! to its derivatives
1745 end subroutine sc_angular_nucl
1747 !-----------------------------------------------------------------------------
1748 subroutine int_bounds(total_ints,lower_bound,upper_bound)
1749 ! implicit real*8 (a-h,o-z)
1750 ! include 'DIMENSIONS'
1752 ! include 'COMMON.SETUP'
1753 integer :: total_ints,lower_bound,upper_bound,nint
1754 integer,dimension(0:nfgtasks) :: int4proc,sint4proc !(0:max_fg_procs)
1755 integer :: i,nexcess
1756 nint=total_ints/nfgtasks
1760 nexcess=total_ints-nint*nfgtasks
1762 int4proc(nfgtasks-i)=int4proc(nfgtasks-i)+1
1766 lower_bound=lower_bound+int4proc(i)
1768 upper_bound=lower_bound+int4proc(fg_rank)
1769 lower_bound=lower_bound+1
1771 end subroutine int_bounds
1772 !-----------------------------------------------------------------------------
1773 subroutine int_bounds1(total_ints,lower_bound,upper_bound)
1774 ! implicit real*8 (a-h,o-z)
1775 ! include 'DIMENSIONS'
1777 ! include 'COMMON.SETUP'
1778 integer :: total_ints,lower_bound,upper_bound,nint
1779 integer :: nexcess,i
1780 integer,dimension(0:nfgtasks) :: int4proc,sint4proc !(0:max_fg_procs)
1781 nint=total_ints/nfgtasks1
1785 nexcess=total_ints-nint*nfgtasks1
1787 int4proc(nfgtasks1-i)=int4proc(nfgtasks1-i)+1
1791 lower_bound=lower_bound+int4proc(i)
1793 upper_bound=lower_bound+int4proc(fg_rank1)
1794 lower_bound=lower_bound+1
1796 end subroutine int_bounds1
1797 !-----------------------------------------------------------------------------
1799 !-----------------------------------------------------------------------------
1800 subroutine chainbuild_cart
1801 ! implicit real*8 (a-h,o-z)
1802 ! include 'DIMENSIONS'
1807 ! include 'COMMON.SETUP'
1808 ! include 'COMMON.CHAIN'
1809 ! include 'COMMON.LOCAL'
1810 ! include 'COMMON.TIME1'
1811 ! include 'COMMON.IOUNITS'
1812 integer :: j,i,ierror,ierr
1813 real(kind=8) :: time00,time01
1815 if (nfgtasks.gt.1) then
1816 ! write (iout,*) "BCAST in chainbuild_cart"
1818 ! Broadcast the order to build the chain and compute internal coordinates
1819 ! to the slaves. The slaves receive the order in ERGASTULUM.
1821 ! write (iout,*) "CHAINBUILD_CART: DC before BCAST"
1823 ! write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1824 ! & (dc(j,i+nres),j=1,3)
1827 call MPI_Bcast(7,1,MPI_INTEGER,king,FG_COMM,IERROR)
1828 time_bcast7=time_bcast7+MPI_Wtime()-time00
1830 call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,&
1832 ! write (iout,*) "CHAINBUILD_CART: DC after BCAST"
1834 ! write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1835 ! & (dc(j,i+nres),j=1,3)
1837 ! write (iout,*) "End BCAST in chainbuild_cart"
1839 time_bcast=time_bcast+MPI_Wtime()-time00
1840 time_bcastc=time_bcastc+MPI_Wtime()-time01
1848 c(j,i)=c(j,i-1)+dc(j,i-1)
1853 c(j,i+nres)=c(j,i)+dc(j,i+nres)
1856 ! write (iout,*) "CHAINBUILD_CART"
1858 call int_from_cart1(.false.)
1860 end subroutine chainbuild_cart
1861 !-----------------------------------------------------------------------------
1863 !-----------------------------------------------------------------------------
1864 real(kind=8) function alpha(i1,i2,i3)
1866 ! Calculates the planar angle between atoms (i1), (i2), and (i3).
1868 ! implicit real*8 (a-h,o-z)
1869 ! include 'DIMENSIONS'
1870 ! include 'COMMON.GEO'
1871 ! include 'COMMON.CHAIN'
1874 real(kind=8) :: x12,x23,y12,y23,z12,z23,vnorm,wnorm,scalar
1881 vnorm=dsqrt(x12*x12+y12*y12+z12*z12)
1882 wnorm=dsqrt(x23*x23+y23*y23+z23*z23)
1883 scalar=(x12*x23+y12*y23+z12*z23)/(vnorm*wnorm)
1887 !-----------------------------------------------------------------------------
1888 real(kind=8) function beta(i1,i2,i3,i4)
1890 ! Calculates the dihedral angle between atoms (i1), (i2), (i3) and (i4)
1892 ! implicit real*8 (a-h,o-z)
1893 ! include 'DIMENSIONS'
1894 ! include 'COMMON.GEO'
1895 ! include 'COMMON.CHAIN'
1897 integer :: i1,i2,i3,i4
1898 real(kind=8) :: x12,x23,x34,y12,y23,y34,z12,z23,z34
1899 real(kind=8) :: wx,wy,wz,wnorm,vx,vy,vz,vnorm,scalar,angle
1900 real(kind=8) :: tx,ty,tz
1910 !d print '(2i3,3f10.5)',i1,i2,x12,y12,z12
1911 !d print '(2i3,3f10.5)',i2,i3,x23,y23,z23
1912 !d print '(2i3,3f10.5)',i3,i4,x34,y34,z34
1916 wnorm=dsqrt(wx*wx+wy*wy+wz*wz)
1920 vnorm=dsqrt(vx*vx+vy*vy+vz*vz)
1921 if (vnorm.gt.1.0D-13 .and. wnorm.gt.1.0D-13) then
1922 scalar=(vx*wx+vy*wy+vz*wz)/(vnorm*wnorm)
1923 if (dabs(scalar).gt.1.0D0) &
1924 scalar=0.99999999999999D0*scalar/dabs(scalar)
1926 !d print '(2i4,10f7.3)',i2,i3,vx,vy,vz,wx,wy,wz,vnorm,wnorm,
1931 ! if (angle.le.0.0D0) angle=pi+angle
1935 scalar=tx*x23+ty*y23+tz*z23
1936 if (scalar.lt.0.0D0) angle=-angle
1940 !-----------------------------------------------------------------------------
1941 real(kind=8) function dist(i1,i2)
1943 ! Calculates the distance between atoms (i1) and (i2).
1945 ! implicit real*8 (a-h,o-z)
1946 ! include 'DIMENSIONS'
1947 ! include 'COMMON.GEO'
1948 ! include 'COMMON.CHAIN'
1951 real(kind=8) :: x12,y12,z12
1955 dist=dsqrt(x12*x12+y12*y12+z12*z12)
1958 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
1959 !-----------------------------------------------------------------------------
1961 !-----------------------------------------------------------------------------
1962 subroutine local_move_init(debug)
1966 ! implicit real*8 (a-h,o-z)
1967 ! include 'DIMENSIONS' ! Needed by COMMON.LOCAL
1968 ! include 'COMMON.GEO' ! For pi, deg2rad
1969 ! include 'COMMON.LOCAL' ! For vbl
1970 ! include 'COMMON.LOCMOVE'
1976 ! Determine wheter to do some debugging output
1977 locmove_output=debug
1979 ! Set the init_called flag to 1
1982 ! The following are never changed
1983 min_theta=60.D0*deg2rad ! (0,PI)
1984 max_theta=175.D0*deg2rad ! (0,PI)
1985 dmin2=vbl*vbl*2.*(1.-cos(min_theta))
1986 dmax2=vbl*vbl*2.*(1.-cos(max_theta))
1989 small2=0.5*small*small
1991 ! Not really necessary...
1997 end subroutine local_move_init
1998 !-----------------------------------------------------------------------------
1999 subroutine local_move(n_start, n_end, PHImin, PHImax)
2000 ! Perform a local move between residues m and n (inclusive)
2001 ! PHImin and PHImax [0,PI] determine the size of the move
2002 ! Works on whatever structure is in the variables theta and phi,
2003 ! sidechain variables are left untouched
2004 ! The final structure is NOT minimized, but both the cartesian
2005 ! variables c and the angles are up-to-date at the end (no further
2006 ! chainbuild is required)
2008 use random,only:ran_number
2010 ! implicit real*8 (a-h,o-z)
2011 ! include 'DIMENSIONS'
2012 ! include 'COMMON.GEO'
2013 ! include 'COMMON.CHAIN'
2014 ! include 'COMMON.VAR'
2015 ! include 'COMMON.MINIM'
2016 ! include 'COMMON.SBRIDGE'
2017 ! include 'COMMON.LOCMOVE'
2019 ! External functions
2020 !EL integer move_res
2021 !EL external move_res
2022 !EL double precision ran_number
2023 !EL external ran_number
2026 integer :: n_start, n_end ! First and last residues to move
2027 real(kind=8) :: PHImin, PHImax ! min/max angles [0,PI]
2031 real(kind=8) :: min,max
2035 ! Check if local_move_init was called. This assumes that it
2036 ! would not be 1 if not explicitely initialized
2037 if (init_called.ne.1) then
2038 write(6,*)' *** local_move_init not called!!!'
2042 ! Quick check for crazy range
2043 if (n_start.gt.n_end .or. n_start.lt.1 .or. n_end.gt.nres) then
2044 write(6,'(a,i3,a,i3)') &
2045 ' *** Cannot make local move between n_start = ',&
2046 n_start,' and n_end = ',n_end
2050 ! Take care of end residues first...
2051 if (n_start.eq.1) then
2052 ! Move residue 1 (completely random)
2053 theta(3)=ran_number(min_theta,max_theta)
2054 phi(4)=ran_number(-PI,PI)
2059 if (n_end.eq.nres) then
2060 ! Move residue nres (completely random)
2061 theta(nres)=ran_number(min_theta,max_theta)
2062 phi(nres)=ran_number(-PI,PI)
2068 ! ...then go through all other residues one by one
2069 ! Start from the two extremes and converge
2074 !$$$c Move the first two residues by less than the others
2075 !$$$ if (i-n_start.lt.3) then
2076 !$$$ if (i-n_start.eq.0) then
2079 !$$$ else if (i-n_start.eq.1) then
2082 !$$$ else if (i-n_start.eq.2) then
2088 ! The actual move, on residue i
2089 iretcode=move_res(min,max,i) ! Discard iretcode
2095 !$$$c Move the last two residues by less than the others
2096 !$$$ if (n_end-j.lt.3) then
2097 !$$$ if (n_end-j.eq.0) then
2100 !$$$ else if (n_end-j.eq.1) then
2103 !$$$ else if (n_end-j.eq.2) then
2109 ! The actual move, on residue j
2110 iretcode=move_res(min,max,j) ! Discard iretcode
2115 call int_from_cart(.false.,.false.)
2118 end subroutine local_move
2119 !-----------------------------------------------------------------------------
2120 subroutine output_tabs
2121 ! Prints out the contents of a_..., b_..., res_...
2125 ! include 'COMMON.GEO'
2126 ! include 'COMMON.LOCMOVE'
2132 write(6,'(8f7.1)')(a_ang(i)*rad2deg,i=0,a_n-1)
2133 write(6,'(8(2x,3l1,2x))')((a_tab(i,j),i=0,2),j=0,a_n-1)
2136 write(6,'(4f7.1)')(b_ang(i)*rad2deg,i=0,b_n-1)
2137 write(6,'(4(2x,3l1,2x))')((b_tab(i,j),i=0,2),j=0,b_n-1)
2140 write(6,'(12f7.1)')(res_ang(i)*rad2deg,i=0,res_n-1)
2141 write(6,'(12(2x,3l1,2x))')((res_tab(0,i,j),i=0,2),j=0,res_n-1)
2142 write(6,'(12(2x,3l1,2x))')((res_tab(1,i,j),i=0,2),j=0,res_n-1)
2143 write(6,'(12(2x,3l1,2x))')((res_tab(2,i,j),i=0,2),j=0,res_n-1)
2146 end subroutine output_tabs
2147 !-----------------------------------------------------------------------------
2148 subroutine angles2tab(PHImin,PHImax,n,ang,tab)
2149 ! Only uses angles if [0,PI] (but PHImin cannot be 0.,
2150 ! and PHImax cannot be PI)
2154 ! include 'COMMON.GEO'
2157 real(kind=8) :: PHImin,PHImax
2161 real(kind=8),dimension(0:3) :: ang
2162 logical,dimension(0:2,0:3) :: tab
2165 if (PHImin .eq. PHImax) then
2166 ! Special case with two 010's
2176 else if (PHImin .eq. PI) then
2177 ! Special case with one 010
2183 else if (PHImax .eq. 0.) then
2184 ! Special case with one 010
2193 if (PHImin .gt. 0.) then
2194 ! Start of range (011)
2199 ! End of range (110)
2203 tab(2,n+1) = .false.
2206 if (PHImax .lt. PI) then
2207 ! Start of range (011)
2212 ! End of range (110)
2216 tab(2,n+1) = .false.
2222 end subroutine angles2tab
2223 !-----------------------------------------------------------------------------
2224 subroutine minmax_angles(x,y,z,r,n,ang,tab)
2225 ! When solutions do not exist, assume all angles
2226 ! are acceptable - i.e., initial geometry must be correct
2230 ! include 'COMMON.GEO'
2231 ! include 'COMMON.LOCMOVE'
2234 real(kind=8) :: x,y,z,r
2238 real(kind=8),dimension(0:3) :: ang
2239 logical,dimension(0:2,0:3) :: tab
2242 real(kind=8) :: num, denom, phi
2243 real(kind=8) :: Kmin, Kmax
2247 num = x*x + y*y + z*z
2250 if (denom .gt. 0.) then
2252 denom = 2.*r*sqrt(denom)
2254 Kmin = (num - dmin2)/denom
2255 Kmax = (num - dmax2)/denom
2257 ! Allowed values of K (else all angles are acceptable)
2260 if (Kmin .gt. 1. .or. abs(Kmin-1.) .lt. small2) then
2262 else if (Kmin .lt. -1. .or. abs(Kmin+1.) .lt. small2) then
2268 if (Kmax .lt. -1. .or. abs(Kmax+1.) .lt. small2) then
2270 else if (Kmax .gt. 1. .or. abs(Kmax-1.) .lt. small2) then
2276 if (Kmax .lt. Kmin) Kmax = Kmin
2278 call angles2tab(Kmin, Kmax, n, ang, tab)
2280 ! Add phi and check that angles are within range (-PI,PI]
2283 if (ang(i) .le. -PI) then
2284 ang(i) = ang(i)+2.*PI
2285 else if (ang(i) .gt. PI) then
2286 ang(i) = ang(i)-2.*PI
2292 end subroutine minmax_angles
2293 !-----------------------------------------------------------------------------
2294 subroutine construct_tab
2295 ! Take a_... and b_... values and produces the results res_...
2296 ! x_ang are assumed to be all different (diff > small)
2297 ! x_tab(1,i) must be 1 for all i (i.e., all x_ang are acceptable)
2301 ! include 'COMMON.LOCMOVE'
2304 integer :: n_max,i,j,index
2310 if (n_max .eq. 0) then
2317 res_tab(j,0,i) = .true.
2318 res_tab(j,2,i) = .true.
2319 res_tab(j,1,i) = .false.
2326 do while (.not.done)
2327 res_ang(index) = flag
2331 if ((a_ang(i)-phi).gt.small .and. &
2332 a_ang(i) .lt. res_ang(index)) then
2333 ! Found a lower angle
2334 res_ang(index) = a_ang(i)
2335 ! Copy the values from a_tab into res_tab(0,,)
2336 res_tab(0,0,index) = a_tab(0,i)
2337 res_tab(0,1,index) = a_tab(1,i)
2338 res_tab(0,2,index) = a_tab(2,i)
2339 ! Set default values for res_tab(1,,)
2340 res_tab(1,0,index) = .true.
2341 res_tab(1,1,index) = .false.
2342 res_tab(1,2,index) = .true.
2343 else if (abs(a_ang(i)-res_ang(index)).lt.small) then
2344 ! Found an equal angle (can only be equal to a b_ang)
2345 res_tab(0,0,index) = a_tab(0,i)
2346 res_tab(0,1,index) = a_tab(1,i)
2347 res_tab(0,2,index) = a_tab(2,i)
2352 if ((b_ang(i)-phi).gt.small .and. &
2353 b_ang(i) .lt. res_ang(index)) then
2354 ! Found a lower angle
2355 res_ang(index) = b_ang(i)
2356 ! Copy the values from b_tab into res_tab(1,,)
2357 res_tab(1,0,index) = b_tab(0,i)
2358 res_tab(1,1,index) = b_tab(1,i)
2359 res_tab(1,2,index) = b_tab(2,i)
2360 ! Set default values for res_tab(0,,)
2361 res_tab(0,0,index) = .true.
2362 res_tab(0,1,index) = .false.
2363 res_tab(0,2,index) = .true.
2364 else if (abs(b_ang(i)-res_ang(index)).lt.small) then
2365 ! Found an equal angle (can only be equal to an a_ang)
2366 res_tab(1,0,index) = b_tab(0,i)
2367 res_tab(1,1,index) = b_tab(1,i)
2368 res_tab(1,2,index) = b_tab(2,i)
2372 if (res_ang(index) .eq. flag) then
2375 else if (index .eq. n_max-1) then
2379 phi = res_ang(index) ! Store previous angle
2387 if (a_n .gt. 0) then
2388 do while (.not.res_tab(0,1,index))
2391 done = res_tab(0,2,index)
2392 do i=index+1,res_n-1
2393 if (res_tab(0,1,i)) then
2394 done = res_tab(0,2,i)
2396 res_tab(0,0,i) = done
2397 res_tab(0,1,i) = done
2398 res_tab(0,2,i) = done
2401 done = res_tab(0,0,index)
2403 if (res_tab(0,1,i)) then
2404 done = res_tab(0,0,i)
2406 res_tab(0,0,i) = done
2407 res_tab(0,1,i) = done
2408 res_tab(0,2,i) = done
2413 res_tab(0,0,i) = .true.
2414 res_tab(0,1,i) = .true.
2415 res_tab(0,2,i) = .true.
2420 if (b_n .gt. 0) then
2421 do while (.not.res_tab(1,1,index))
2424 done = res_tab(1,2,index)
2425 do i=index+1,res_n-1
2426 if (res_tab(1,1,i)) then
2427 done = res_tab(1,2,i)
2429 res_tab(1,0,i) = done
2430 res_tab(1,1,i) = done
2431 res_tab(1,2,i) = done
2434 done = res_tab(1,0,index)
2436 if (res_tab(1,1,i)) then
2437 done = res_tab(1,0,i)
2439 res_tab(1,0,i) = done
2440 res_tab(1,1,i) = done
2441 res_tab(1,2,i) = done
2446 res_tab(1,0,i) = .true.
2447 res_tab(1,1,i) = .true.
2448 res_tab(1,2,i) = .true.
2452 ! Finally fill the last row with AND operation
2455 res_tab(2,j,i) = (res_tab(0,j,i) .and. res_tab(1,j,i))
2460 end subroutine construct_tab
2461 !-----------------------------------------------------------------------------
2462 subroutine construct_ranges(phi_n,phi_start,phi_end)
2463 ! Given the data in res_..., construct a table of
2464 ! min/max allowed angles
2468 ! include 'COMMON.GEO'
2469 ! include 'COMMON.LOCMOVE'
2473 real(kind=8),dimension(0:11) :: phi_start,phi_end
2480 if (res_n .eq. 0) then
2481 ! Any move is allowed
2489 do while (.not.done)
2490 ! Find start of range (01x)
2492 do while (.not.done)
2493 if (res_tab(2,0,index).or.(.not.res_tab(2,1,index))) then
2497 phi_start(phi_n) = res_ang(index)
2499 if (index .eq. res_n) done = .true.
2501 ! If a start was found (index < res_n), find the end of range (x10)
2502 ! It may not be found without wrapping around
2503 if (index .lt. res_n) then
2505 do while (.not.done)
2506 if ((.not.res_tab(2,1,index)).or.res_tab(2,2,index)) then
2511 if (index .eq. res_n) done = .true.
2513 if (index .lt. res_n) then
2514 ! Found the end of the range
2515 phi_end(phi_n) = res_ang(index)
2518 if (index .eq. res_n) then
2524 ! Need to wrap around
2526 phi_end(phi_n) = flag
2530 ! Take care of the last one if need to wrap around
2531 if (phi_end(phi_n) .eq. flag) then
2533 do while ((.not.res_tab(2,1,index)).or.res_tab(2,2,index))
2536 phi_end(phi_n) = res_ang(index) + 2.*PI
2542 end subroutine construct_ranges
2543 !-----------------------------------------------------------------------------
2544 subroutine fix_no_moves(phi)
2548 ! include 'COMMON.GEO'
2549 ! include 'COMMON.LOCMOVE'
2556 real(kind=8) :: diff,temp
2559 ! Look for first 01x in gammas (there MUST be at least one)
2562 do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index)))
2565 if (res_ang(index) .le. 0.D0) then ! Make sure it's from PHImax
2566 ! Try to increase PHImax
2567 if (index .gt. 0) then
2568 phi = res_ang(index-1)
2569 diff = abs(res_ang(index) - res_ang(index-1))
2571 ! Look for last (corresponding) x10
2573 do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index))
2576 if (index .lt. res_n-1) then
2577 temp = abs(res_ang(index) - res_ang(index+1))
2578 if (temp .lt. diff) then
2579 phi = res_ang(index+1)
2585 ! If increasing PHImax didn't work, decreasing PHImin
2586 ! will (with one exception)
2587 ! Look for first x10 (there MUST be at least one)
2589 do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index))
2592 if (res_ang(index) .lt. 0.D0) then ! Make sure it's from PHImin
2593 ! Try to decrease PHImin
2594 if (index .lt. res_n-1) then
2595 temp = abs(res_ang(index) - res_ang(index+1))
2596 if (res_ang(index+1) .le. 0.D0 .and. temp .lt. diff) then
2597 phi = res_ang(index+1)
2601 ! Look for last (corresponding) 01x
2603 do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index)))
2606 if (index .gt. 0) then
2607 temp = abs(res_ang(index) - res_ang(index-1))
2608 if (res_ang(index-1) .ge. 0.D0 .and. temp .lt. diff) then
2609 phi = res_ang(index-1)
2615 ! If it still didn't work, it must be PHImax == 0. or PHImin == PI
2616 if (diff .eq. flag) then
2618 if (res_tab(index,1,0) .or. (.not.res_tab(index,1,1)) .or. &
2619 res_tab(index,1,2)) index = res_n - 1
2620 ! This MUST work at this point
2621 if (index .eq. 0) then
2624 phi = res_ang(index - 1)
2629 end subroutine fix_no_moves
2630 !-----------------------------------------------------------------------------
2631 integer function move_res(PHImin,PHImax,i_move)
2632 ! Moves residue i_move (in array c), leaving everything else fixed
2633 ! Starting geometry is not checked, it should be correct!
2634 ! R(,i_move) is the only residue that will move, but must have
2635 ! 1 < i_move < nres (i.e., cannot move ends)
2636 ! Whether any output is done is controlled by locmove_output
2638 use random,only:ran_number
2640 ! implicit real*8 (a-h,o-z)
2641 ! include 'DIMENSIONS'
2642 ! include 'COMMON.CHAIN'
2643 ! include 'COMMON.GEO'
2644 ! include 'COMMON.LOCMOVE'
2646 ! External functions
2647 !EL double precision ran_number
2648 !EL external ran_number
2651 real(kind=8) :: PHImin,PHImax
2655 ! 0: move successfull
2656 ! 1: Dmin or Dmax had to be modified
2657 ! 2: move failed - check your input geometry
2661 real(kind=8),dimension(0:2) :: X,Y,Z,Orig
2662 real(kind=8),dimension(0:2) :: P
2663 logical :: no_moves,done
2664 integer :: index,i,j
2665 real(kind=8) :: phi,temp,radius
2666 real(kind=8),dimension(0:11) :: phi_start,phi_end
2669 ! Set up the coordinate system
2671 Orig(i)=0.5*(c(i+1,i_move-1)+c(i+1,i_move+1)) ! Position of origin
2675 Z(i)=c(i+1,i_move+1)-c(i+1,i_move-1)
2677 temp=sqrt(Z(0)*Z(0)+Z(1)*Z(1)+Z(2)*Z(2))
2683 X(i)=c(i+1,i_move)-Orig(i)
2685 ! radius is the radius of the circle on which c(,i_move) can move
2686 radius=sqrt(X(0)*X(0)+X(1)*X(1)+X(2)*X(2))
2691 Y(0)=Z(1)*X(2)-X(1)*Z(2)
2692 Y(1)=X(0)*Z(2)-Z(0)*X(2)
2693 Y(2)=Z(0)*X(1)-X(0)*Z(1)
2695 ! Calculate min, max angles coming from dmin, dmax to c(,i_move-2)
2696 if (i_move.gt.2) then
2698 P(i)=c(i+1,i_move-2)-Orig(i)
2700 call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2),&
2701 P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2),&
2702 P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2),&
2703 radius,a_n,a_ang,a_tab)
2708 ! Calculate min, max angles coming from dmin, dmax to c(,i_move+2)
2709 if (i_move.lt.nres-2) then
2711 P(i)=c(i+1,i_move+2)-Orig(i)
2713 call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2),&
2714 P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2),&
2715 P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2),&
2716 radius,b_n,b_ang,b_tab)
2721 ! Construct the resulting table for alpha and beta
2722 call construct_tab()
2724 if (locmove_output) then
2725 print *,'ALPHAS & BETAS TABLE'
2729 ! Check that there is at least one possible move
2731 if (res_n .eq. 0) then
2735 do while ((index .lt. res_n) .and. no_moves)
2736 if (res_tab(2,1,index)) no_moves = .false.
2741 if (locmove_output) print *,' *** Cannot move anywhere'
2746 ! Transfer res_... into a_...
2749 if ( (res_tab(2,0,i).neqv.res_tab(2,1,i)) .or. &
2750 (res_tab(2,0,i).neqv.res_tab(2,2,i)) ) then
2751 a_ang(a_n) = res_ang(i)
2753 a_tab(j,a_n) = res_tab(2,j,i)
2759 ! Check that the PHI's are within [0,PI]
2760 if (PHImin .lt. 0. .or. abs(PHImin) .lt. small) PHImin = -flag
2761 if (PHImin .gt. PI .or. abs(PHImin-PI) .lt. small) PHImin = PI
2762 if (PHImax .gt. PI .or. abs(PHImax-PI) .lt. small) PHImax = flag
2763 if (PHImax .lt. 0. .or. abs(PHImax) .lt. small) PHImax = 0.
2764 if (PHImax .lt. PHImin) PHImax = PHImin
2765 ! Calculate min and max angles coming from PHImin and PHImax,
2766 ! and put them in b_...
2767 call angles2tab(PHImin, PHImax, b_n, b_ang, b_tab)
2768 ! Construct the final table
2769 call construct_tab()
2771 if (locmove_output) then
2772 print *,'FINAL TABLE'
2776 ! Check that there is at least one possible move
2778 if (res_n .eq. 0) then
2782 do while ((index .lt. res_n) .and. no_moves)
2783 if (res_tab(2,1,index)) no_moves = .false.
2789 ! Take care of the case where no solution exists...
2790 call fix_no_moves(phi)
2791 if (locmove_output) then
2792 print *,' *** Had to modify PHImin or PHImax'
2793 print *,'phi: ',phi*rad2deg
2797 ! ...or calculate the solution
2798 ! Construct phi_start/phi_end arrays
2799 call construct_ranges(phi_n, phi_start, phi_end)
2800 ! Choose random angle phi in allowed range(s)
2803 temp = temp + phi_end(i) - phi_start(i)
2805 phi = ran_number(phi_start(0),phi_start(0)+temp)
2808 do while (.not.done)
2809 if (phi .lt. phi_end(index)) then
2814 if (index .eq. phi_n) then
2816 else if (.not.done) then
2817 phi = phi + phi_start(index) - phi_end(index-1)
2820 if (index.eq.phi_n) phi=phi_end(phi_n-1) ! Fix numerical errors
2821 if (phi .gt. PI) phi = phi-2.*PI
2823 if (locmove_output) then
2824 print *,'ALLOWED RANGE(S)'
2826 print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg
2828 print *,'phi: ',phi*rad2deg
2833 ! Re-use radius as temp variable
2834 temp=radius*cos(phi)
2835 radius=radius*sin(phi)
2837 c(i+1,i_move)=Orig(i)+temp*X(i)+radius*Y(i)
2841 end function move_res
2842 !-----------------------------------------------------------------------------
2847 ! implicit real*8 (a-h,o-z)
2848 ! include 'DIMENSIONS'
2849 ! include 'COMMON.GEO'
2850 ! include 'COMMON.LOCAL'
2851 ! include 'COMMON.LOCMOVE'
2853 ! External functions
2854 !EL integer move_res
2855 !EL external move_res
2860 real(kind=8),dimension(0:11) :: phi_start,phi_end
2862 real(kind=8),dimension(0:2,0:5) :: R
2864 locmove_output=.true.
2866 ! call angles2tab(30.*deg2rad,70.*deg2rad,a_n,a_ang,a_tab)
2867 ! call angles2tab(80.*deg2rad,130.*deg2rad,b_n,b_ang,b_tab)
2868 ! call minmax_angles(0.D0,3.8D0,0.D0,3.8D0,b_n,b_ang,b_tab)
2869 ! call construct_tab
2872 ! call construct_ranges(phi_n,phi_start,phi_end)
2874 ! print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg
2877 ! call fix_no_moves(phi)
2878 ! print *,'NO MOVES FOUND, BEST PHI IS',phi*rad2deg
2884 R(1,1)=-cos(28.D0*deg2rad)
2885 R(2,1)=-0.5D0-sin(28.D0*deg2rad)
2889 R(0,3)=cos(30.D0*deg2rad)
2896 R(1,5)=cos(26.D0*deg2rad)
2897 R(2,5)=0.5D0+sin(26.D0*deg2rad)
2903 ! i=move_res(R(0,1),0.D0*deg2rad,180.D0*deg2rad)
2905 i=move_res(0.D0*deg2rad,180.D0*deg2rad,imov)
2906 print *,'RETURNED ',i
2907 print *,(R(i,3)/vbl,i=0,2)
2910 end subroutine loc_test
2912 !-----------------------------------------------------------------------------
2914 !-----------------------------------------------------------------------------
2915 subroutine MATMULT(A1,A2,A3)
2916 ! implicit real*8 (a-h,o-z)
2917 ! include 'DIMENSIONS'
2920 real(kind=8) :: A3IJ
2922 real(kind=8),DIMENSION(3,3) :: A1,A2,A3
2923 real(kind=8),DIMENSION(3,3) :: AI3
2928 3 A3IJ=A3IJ+A1(I,K)*A2(K,J)
2936 end subroutine MATMULT
2937 !-----------------------------------------------------------------------------
2939 !-----------------------------------------------------------------------------
2940 subroutine int_from_cart(lside,lprn)
2941 ! implicit real*8 (a-h,o-z)
2942 ! include 'DIMENSIONS'
2943 use control_data,only:out1file
2947 ! include 'COMMON.LOCAL'
2948 ! include 'COMMON.VAR'
2949 ! include 'COMMON.CHAIN'
2950 ! include 'COMMON.INTERACT'
2951 ! include 'COMMON.IOUNITS'
2952 ! include 'COMMON.GEO'
2953 ! include 'COMMON.NAMES'
2954 ! include 'COMMON.CONTROL'
2955 ! include 'COMMON.SETUP'
2956 character(len=3) :: seq,res
2958 character(len=80) :: card
2959 real(kind=8),dimension(3,20) :: sccor
2960 integer :: i,j,iti !el rescode,
2961 logical :: lside,lprn
2962 real(kind=8) :: di,cosfac,sinfac
2966 if(me.eq.king.or..not.out1file)then
2968 write (iout,'(/a)') &
2969 'Internal coordinates calculated from crystal structure.'
2971 write (iout,'(8a)') ' Res ',' dvb',' Theta',&
2972 ' Gamma',' Dsc_id',' Dsc',' Alpha',&
2975 write (iout,'(4a)') ' Res ',' dvb',' Theta',&
2981 ! if (molnum(i).ne.1) cycle
2982 !in wham do i=1,nres
2984 if (((dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0).and.&
2985 (iti.ne.ntyp1 .and. itype(i+1,1).ne.ntyp1)).and.molnum(i).eq.1) then
2986 write (iout,'(a,i4)') 'Bad Cartesians for residue',i
2990 vbld(i+1)=dist(i,i+1)
2991 vbld_inv(i+1)=1.0d0/vbld(i+1)
2993 if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
2994 if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
2998 ! if (itype(1,1).eq.ntyp1) then
3000 ! c(j,1)=c(j,2)+(c(j,3)-c(j,4))
3003 ! if (itype(nres,1).eq.ntyp1) then
3005 ! c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
3009 ! if (unres_pdb) then
3010 ! if (itype(1,1).eq.21) then
3011 ! theta(3)=90.0d0*deg2rad
3012 ! phi(4)=180.0d0*deg2rad
3014 ! vbld_inv(2)=1.0d0/vbld(2)
3016 ! if (itype(nres,1).eq.21) then
3017 ! theta(nres)=90.0d0*deg2rad
3018 ! phi(nres)=180.0d0*deg2rad
3020 ! vbld_inv(nres)=1.0d0/vbld(2)
3026 c(j,nres2+2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i) &
3027 +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
3028 ! in wham c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1)
3033 ! 10/03/12 Adam: Correction for zero SC-SC bond length
3035 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1 .and. di.eq.0.0d0) &
3036 di=dsc(itype(i,molnum(i)))
3038 if (itype(i,1).ne.10) then
3039 vbld_inv(i+nres)=1.0d0/di
3041 vbld_inv(i+nres)=0.0d0
3045 alph(i)=alpha(nres+i,i,nres2+2)
3046 omeg(i)=beta(nres+i,i,nres2+2,i+1)
3049 if(me.eq.king.or..not.out1file)then
3051 write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,vbld(i),&
3052 rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),&
3053 rad2deg*alph(i),rad2deg*omeg(i)
3056 if(me.eq.king.or..not.out1file)then
3058 write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,vbld(i),&
3059 rad2deg*theta(i),rad2deg*phi(i),dsc(iti+1),vbld(nres+i),&
3060 rad2deg*alph(i),rad2deg*omeg(i)
3067 if(me.eq.king.or..not.out1file) &
3068 write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,dist(i,i-1),&
3069 rad2deg*theta(i),rad2deg*phi(i)
3073 end subroutine int_from_cart
3074 !-----------------------------------------------------------------------------
3075 subroutine sc_loc_geom(lprn)
3076 ! implicit real*8 (a-h,o-z)
3077 ! include 'DIMENSIONS'
3078 use control_data,only:out1file
3082 ! include 'COMMON.LOCAL'
3083 ! include 'COMMON.VAR'
3084 ! include 'COMMON.CHAIN'
3085 ! include 'COMMON.INTERACT'
3086 ! include 'COMMON.IOUNITS'
3087 ! include 'COMMON.GEO'
3088 ! include 'COMMON.NAMES'
3089 ! include 'COMMON.CONTROL'
3090 ! include 'COMMON.SETUP'
3091 real(kind=8),dimension(3) :: x_prime,y_prime,z_prime
3094 integer :: i,j,it,iti
3095 real(kind=8) :: cosfac2,sinfac2,xx,yy,zz,cosfac,sinfac
3098 dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
3102 if (itype(i,1).ne.10) then
3104 dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
3108 dc_norm(j,i+nres)=0.0d0
3113 costtab(i+1) =dcos(theta(i+1))
3114 sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
3115 cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
3116 sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
3117 cosfac2=0.5d0/(1.0d0+costtab(i+1))
3118 cosfac=dsqrt(cosfac2)
3119 sinfac2=0.5d0/(1.0d0-costtab(i+1))
3120 sinfac=dsqrt(sinfac2)
3123 if ((it.ne.10).and.(it.ne.ntyp1)) then
3124 !el if (it.ne.10) then
3126 ! Compute the axes of tghe local cartesian coordinates system; store in
3127 ! x_prime, y_prime and z_prime
3135 x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
3136 y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
3138 call vecpr(x_prime,y_prime,z_prime)
3140 ! Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
3141 ! to local coordinate system. Store in xx, yy, zz.
3147 xx = xx + x_prime(j)*dc_norm(j,i+nres)
3148 yy = yy + y_prime(j)*dc_norm(j,i+nres)
3149 zz = zz + z_prime(j)*dc_norm(j,i+nres)
3164 if(me.eq.king.or..not.out1file) &
3165 write (iout,'(a3,i4,3f10.5)') restyp(iti,1),i,xxref(i),&
3171 end subroutine sc_loc_geom
3172 !-----------------------------------------------------------------------------
3173 subroutine sccenter(ires,nscat,sccor)
3174 ! implicit real*8 (a-h,o-z)
3175 ! include 'DIMENSIONS'
3176 ! include 'COMMON.CHAIN'
3177 integer :: i,j,ires,nscat
3178 real(kind=8),dimension(3,20) :: sccor
3179 real(kind=8) :: sccmj
3180 ! print *,"I am in sccenter",ires,nscat
3184 sccmj=sccmj+sccor(j,i)
3185 !C print *,"insccent", ires,sccor(j,i)
3187 dc(j,ires)=sccmj/nscat
3190 end subroutine sccenter
3191 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3192 !-----------------------------------------------------------------------------
3193 subroutine bond_regular
3195 ! implicit real*8 (a-h,o-z)
3196 ! include 'DIMENSIONS'
3197 ! include 'COMMON.VAR'
3198 ! include 'COMMON.LOCAL'
3199 ! include 'COMMON.CALC'
3200 ! include 'COMMON.INTERACT'
3201 ! include 'COMMON.CHAIN'
3205 vbld_inv(i+1)=1.0d0/vbld(i+1)
3206 vbld(i+1+nres)=dsc(itype(i+1,molnum(i)))
3207 vbld_inv(i+1+nres)=dsc_inv(itype(i+1,molnum(i)))
3208 ! print *,vbld(i+1),vbld(i+1+nres)
3211 end subroutine bond_regular
3213 !-----------------------------------------------------------------------------
3215 !-----------------------------------------------------------------------------
3216 subroutine refsys(i2,i3,i4,e1,e2,e3,fail)
3217 ! This subroutine calculates unit vectors of a local reference system
3218 ! defined by atoms (i2), (i3), and (i4). The x axis is the axis from
3219 ! atom (i3) to atom (i2), and the xy plane is the plane defined by atoms
3220 ! (i2), (i3), and (i4). z axis is directed according to the sign of the
3221 ! vector product (i3)-(i2) and (i3)-(i4). Sets fail to .true. if atoms
3222 ! (i2) and (i3) or (i3) and (i4) coincide or atoms (i2), (i3), and (i4)
3223 ! form a linear fragment. Returns vectors e1, e2, and e3.
3224 ! implicit real*8 (a-h,o-z)
3225 ! include 'DIMENSIONS'
3227 real(kind=8),dimension(3) :: e1,e2,e3
3228 real(kind=8),dimension(3) :: u,z
3229 ! include 'COMMON.IOUNITS'
3230 ! include 'COMMON.CHAIN'
3231 real(kind=8) :: coinc=1.0D-13,align=1.0D-13
3233 integer :: i,i1,i2,i3,i4
3234 real(kind=8) :: v1,v2,v3,s1,s2,zi,ui,anorm
3247 if (s1.gt.coinc) goto 2
3248 write (iout,1000) i2,i3,i1
3253 2 if (s2.gt.coinc) goto 4
3254 write(iout,1000) i3,i4,i1
3261 v1=z(2)*u(3)-z(3)*u(2)
3262 v2=z(3)*u(1)-z(1)*u(3)
3263 v3=z(1)*u(2)-z(2)*u(1)
3264 anorm=dsqrt(v1*v1+v2*v2+v3*v3)
3265 if (anorm.gt.align) goto 6
3266 write (iout,1010) i2,i3,i4,i1
3278 e2(1)=e1(3)*e3(2)-e1(2)*e3(3)
3279 e2(2)=e1(1)*e3(3)-e1(3)*e3(1)
3280 e2(3)=e1(2)*e3(1)-e1(1)*e3(2)
3281 1000 format (/1x,' * * * Error - atoms',i4,' and',i4,' coincide.',&
3282 'coordinates of atom',i4,' are set to zero.')
3283 1010 format (/1x,' * * * Error - atoms',2(i4,2h, ),i4,' form a linear',&
3284 ' fragment. coordinates of atom',i4,' are set to zero.')
3286 end subroutine refsys
3287 !-----------------------------------------------------------------------------
3289 !-----------------------------------------------------------------------------
3290 subroutine int_to_cart
3291 !--------------------------------------------------------------
3292 ! This subroutine converts the energy derivatives from internal
3293 ! coordinates to cartesian coordinates
3294 !-------------------------------------------------------------
3295 ! implicit real*8 (a-h,o-z)
3296 ! include 'DIMENSIONS'
3297 ! include 'COMMON.VAR'
3298 ! include 'COMMON.CHAIN'
3299 ! include 'COMMON.DERIV'
3300 ! include 'COMMON.GEO'
3301 ! include 'COMMON.LOCAL'
3302 ! include 'COMMON.INTERACT'
3303 ! include 'COMMON.MD'
3304 ! include 'COMMON.IOUNITS'
3305 ! include 'COMMON.SCCOR'
3306 ! calculating dE/ddc1
3309 ! print *,"gloc",gloc(:,:)
3310 ! print *, "gcart",gcart(:,:)
3311 if (nres.lt.3) go to 18
3313 gcart(j,1)=gcart(j,1)+gloc(1,icg)*dphi(j,1,4) &
3314 +gloc(nres-2,icg)*dtheta(j,1,3)
3315 ! write(iout,*) "pierwszy gcart", gcart(j,2)
3316 if ((itype(2,1).ne.10).and.&
3317 (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)).and.(molnum(2).ne.5))) then
3318 gcart(j,1)=gcart(j,1)+gloc(ialph(2,1),icg)*dalpha(j,1,2)+ &
3319 gloc(ialph(2,1)+nside,icg)*domega(j,1,2)
3322 ! Calculating the remainder of dE/ddc2
3324 gcart(j,2)=gcart(j,2)+gloc(1,icg)*dphi(j,2,4)+ &
3325 gloc(nres-2,icg)*dtheta(j,2,3)+gloc(nres-1,icg)*dtheta(j,1,4)
3326 ! write(iout,*) "drugi gcart", gcart(j,2)
3327 if((itype(2,1).ne.10).and.(molnum(2).ne.5)) then
3328 gcart(j,2)=gcart(j,2)+gloc(ialph(2,1),icg)*dalpha(j,2,2)+ &
3329 gloc(ialph(2,1)+nside,icg)*domega(j,2,2)
3331 if(itype(3,1).ne.10) then
3332 gcart(j,2)=gcart(j,2)+gloc(ialph(3,1),icg)*dalpha(j,1,3)+ &
3333 gloc(ialph(3,1)+nside,icg)*domega(j,1,3)
3336 gcart(j,2)=gcart(j,2)+gloc(2,icg)*dphi(j,1,5)
3339 ! If there are only five residues
3342 gcart(j,3)=gcart(j,3)+gloc(1,icg)*dphi(j,3,4)+gloc(2,icg)* &
3343 dphi(j,2,5)+gloc(nres-1,icg)*dtheta(j,2,4)+gloc(nres,icg)* &
3345 ! if(itype(3,1).ne.10) then
3346 if ((itype(3,1).ne.10).and.&
3347 (itype(3,molnum(3)).ne.ntyp1_molec(molnum(3)))) then
3348 gcart(j,3)=gcart(j,3)+gloc(ialph(3,1),icg)* &
3349 dalpha(j,2,3)+gloc(ialph(3,1)+nside,icg)*domega(j,2,3)
3351 ! if(itype(4,1).ne.10) then
3352 if ((itype(4,1).ne.10).and.&
3353 (itype(4,molnum(4)).ne.ntyp1_molec(molnum(4)))) then
3354 gcart(j,3)=gcart(j,3)+gloc(ialph(4,1),icg)* &
3355 dalpha(j,1,4)+gloc(ialph(4,1)+nside,icg)*domega(j,1,4)
3359 ! If there are more than five residues
3363 gcart(j,i)=gcart(j,i)+gloc(i-2,icg)*dphi(j,3,i+1) &
3364 +gloc(i-1,icg)*dphi(j,2,i+2)+ &
3365 gloc(i,icg)*dphi(j,1,i+3)+gloc(nres+i-4,icg)*dtheta(j,2,i+1)+ &
3366 gloc(nres+i-3,icg)*dtheta(j,1,i+2)
3367 if(itype(i,1).ne.10) then
3368 gcart(j,i)=gcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,2,i)+ &
3369 gloc(ialph(i,1)+nside,icg)*domega(j,2,i)
3371 if(itype(i+1,1).ne.10) then
3372 gcart(j,i)=gcart(j,i)+gloc(ialph(i+1,1),icg)*dalpha(j,1,i+1) &
3373 +gloc(ialph(i+1,1)+nside,icg)*domega(j,1,i+1)
3378 ! Setting dE/ddnres-2
3381 gcart(j,nres-2)=gcart(j,nres-2)+gloc(nres-4,icg)* &
3382 dphi(j,3,nres-1)+gloc(nres-3,icg)*dphi(j,2,nres) &
3383 +gloc(2*nres-6,icg)* &
3384 dtheta(j,2,nres-1)+gloc(2*nres-5,icg)*dtheta(j,1,nres)
3385 if(itype(nres-2,1).ne.10) then
3386 gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-2,1),icg)* &
3387 dalpha(j,2,nres-2)+gloc(ialph(nres-2,1)+nside,icg)* &
3390 if(itype(nres-1,1).ne.10) then
3391 gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-1,1),icg)* &
3392 dalpha(j,1,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
3397 ! Settind dE/ddnres-1
3401 write(iout,*)"in int to carta",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), &
3402 gloc(2*nres-5,icg),dtheta(j,2,nres)
3408 gcart(j,nres-1)=gcart(j,nres-1)+gloc(nres-3,icg)*dphi(j,3,nres)+ &
3409 gloc(2*nres-5,icg)*dtheta(j,2,nres)
3412 write(iout,*)"in int to cartb",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), &
3413 gloc(2*nres-5,icg),dtheta(j,2,nres)
3417 if(itype(nres-1,1).ne.10) then
3418 gcart(j,nres-1)=gcart(j,nres-1)+gloc(ialph(nres-1,1),icg)* &
3419 dalpha(j,2,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
3423 write(iout,*)"in int to cart2",i,gcart(j,nres-1),gloc(ialph(nres-1,1),icg)* &
3424 dalpha(j,2,nres-1),gloc(ialph(nres-1,1)+nside,icg), &
3432 ! The side-chain vector derivatives
3434 if(itype(i,1).ne.10 .and. &
3435 itype(i,molnum(i)).ne.ntyp1_molec(molnum(i)).and.&
3436 (molnum(i).ne.5)) then
3438 gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i) &
3439 +gloc(ialph(i,1)+nside,icg)*domega(j,3,i)
3442 write(iout,*)"in int to cart",i, gxcart(j,i),gloc(ialph(i,1),icg),dalpha(j,3,i), &
3443 gloc(ialph(i,1)+nside,icg),domega(j,3,i)
3449 !----------------------------------------------------------------------
3450 ! INTERTYP=1 SC...Ca...Ca...Ca
3451 ! INTERTYP=2 Ca...Ca...Ca...SC
3452 ! INTERTYP=3 SC...Ca...Ca...SC
3453 ! calculating dE/ddc1
3455 ! write(iout,*) "przed sccor gcart", gcart(1,2)
3459 ! write (iout,*) "poczotkoawy",i,gloc_sc(1,i,icg)
3461 if (nres.lt.2) return
3462 if ((nres.lt.3).and.(itype(1,1).eq.10)) return
3463 if ((itype(1,1).ne.10).and. &
3464 (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))).and.(molnum(1).ne.5)) then
3466 !c Derviative was calculated for oposite vector of side chain therefore
3467 ! there is "-" sign before gloc_sc
3468 gxcart(j,1)=gxcart(j,1)-gloc_sc(1,0,icg)* &
3470 gcart(j,1)=gcart(j,1)+gloc_sc(1,0,icg)* &
3472 if ((itype(2,1).ne.10).and. &
3473 (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2))).and.(molnum(2).ne.5)) then
3474 gxcart(j,1)= gxcart(j,1) &
3475 -gloc_sc(3,0,icg)*dtauangle(j,3,1,3)
3476 gcart(j,1)=gcart(j,1)+gloc_sc(3,0,icg)* &
3481 if ((nres.ge.3).and.(itype(3,molnum(3)).ne.10).and.&
3482 (itype(3,molnum(3)).ne.ntyp1_molec(molnum(3)))) &
3485 gcart(j,1)=gcart(j,1)+gloc_sc(2,1,icg)*dtauangle(j,2,1,4)
3488 ! As potetnial DO NOT depend on omicron anlge their derivative is
3490 ! & +gloc_sc(intertyp,nres-2,icg)*dtheta(j,1,3)
3492 ! write(iout,*) "przed dE/ddc2 gcart", gcart(1,2)
3494 ! Calculating the remainder of dE/ddc2
3496 if((itype(2,1).ne.10).and. &
3497 (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)).and.(molnum(2).ne.5))) then
3498 if ((itype(1,1).ne.10).and.&
3499 ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))).and.(molnum(1).ne.5))&
3500 gxcart(j,2)=gxcart(j,2)+ &
3501 gloc_sc(3,0,icg)*dtauangle(j,3,3,3)
3502 if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,molnum(3)).ne.ntyp1_molec(3)).and.molnum(3).ne.5) &
3504 gxcart(j,2)=gxcart(j,2)-gloc_sc(3,1,icg)*dtauangle(j,3,1,4)
3505 !c the - above is due to different vector direction
3506 gcart(j,2)=gcart(j,2)+gloc_sc(3,1,icg)*dtauangle(j,3,2,4)
3509 ! if ((itype(1,1).ne.10).and.&
3510 ! ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))))) &
3511 gxcart(j,2)=gxcart(j,2)-gloc_sc(1,1,icg)*dtauangle(j,1,1,4)
3512 !c the - above is due to different vector direction
3513 gcart(j,2)=gcart(j,2)+gloc_sc(1,1,icg)*dtauangle(j,1,2,4)
3514 ! write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,2,4),"gcart",gcart(j,2)
3515 ! write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx"
3518 if ((itype(1,1).ne.10).and.&
3519 (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))).and.(molnum(1).ne.5)) then
3520 gcart(j,2)=gcart(j,2)+gloc_sc(1,0,icg)*dtauangle(j,1,3,3)
3521 ! write(iout,*) gloc_sc(1,0,icg),dtauangle(j,1,3,3)
3523 if ((itype(3,1).ne.10).and.(nres.ge.3).and.(molnum(3).ne.5)) then
3524 gcart(j,2)=gcart(j,2)+gloc_sc(2,1,icg)*dtauangle(j,2,2,4)
3525 ! write(iout,*) gloc_sc(2,1,icg),dtauangle(j,2,2,4)
3527 if ((itype(4,1).ne.10).and.(nres.ge.4).and.(molnum(4).ne.5)) then
3528 gcart(j,2)=gcart(j,2)+gloc_sc(2,2,icg)*dtauangle(j,2,1,5)
3529 ! write(iout,*) gloc_sc(2,2,icg),dtauangle(j,2,1,5)
3532 ! write(iout,*) gcart(j,2),itype(2,1),itype(1,1),itype(3,1), "gcart2"
3534 ! write(iout,*) "po dE/ddc2 gcart", gcart(1,2)
3536 ! If there are more than five residues
3540 ! write(iout,*) "before", gcart(j,i)
3541 if ((itype(i,1).ne.10).and.&
3542 (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))).and.(molnum(i).ne.5)) then
3543 gxcart(j,i)=gxcart(j,i)+gloc_sc(2,i-2,icg) &
3544 *dtauangle(j,2,3,i+1) &
3545 -gloc_sc(1,i-1,icg)*dtauangle(j,1,1,i+2)
3546 gcart(j,i)=gcart(j,i)+gloc_sc(1,i-1,icg) &
3547 *dtauangle(j,1,2,i+2)
3548 ! write(iout,*) "new",j,i,
3549 ! & gcart(j,i),gloc_sc(1,i-1,icg),dtauangle(j,1,2,i+2)
3550 ! if (itype(i-1,1).ne.10) then
3551 if ((itype(i-1,1).ne.10).and.&
3552 (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1))).and.(molnum(i-1).eq.5)) then
3554 gxcart(j,i)=gxcart(j,i)+gloc_sc(3,i-2,icg) &
3555 *dtauangle(j,3,3,i+1)
3557 ! if (itype(i+1,1).ne.10) then
3558 if ((itype(i+1,1).ne.10).and.&
3559 (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1)))) then
3560 gxcart(j,i)=gxcart(j,i)-gloc_sc(3,i-1,icg) &
3561 *dtauangle(j,3,1,i+2)
3562 gcart(j,i)=gcart(j,i)+gloc_sc(3,i-1,icg) &
3563 *dtauangle(j,3,2,i+2)
3566 ! if (itype(i-1,1).ne.10) then
3567 if ((itype(i-1,1).ne.10).and.&
3568 (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1)))) then
3569 gcart(j,i)=gcart(j,i)+gloc_sc(1,i-2,icg)* &
3570 dtauangle(j,1,3,i+1)
3572 ! if (itype(i+1,1).ne.10) then
3573 if ((itype(i+1,1).ne.10).and.&
3574 (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1)))) then
3575 gcart(j,i)=gcart(j,i)+gloc_sc(2,i-1,icg)* &
3576 dtauangle(j,2,2,i+2)
3577 ! write(iout,*) "numer",i,gloc_sc(2,i-1,icg),
3578 ! & dtauangle(j,2,2,i+2)
3580 ! if (itype(i+2,1).ne.10) then
3581 if ((itype(i+2,1).ne.10).and.&
3582 (itype(i+2,molnum(i+2)).ne.ntyp1_molec(molnum(i+2))).and.(molnum(i+2).ne.5)) then
3583 gcart(j,i)=gcart(j,i)+gloc_sc(2,i,icg)* &
3584 dtauangle(j,2,1,i+3)
3589 ! Setting dE/ddnres-1
3592 if ((itype(nres-1,1).ne.10).and.&
3593 (itype(nres-1,molnum(nres-1)).ne.ntyp1_molec(molnum(nres-1))).and.(molnum(nres-1).ne.5)) then
3594 gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(2,nres-3,icg) &
3595 *dtauangle(j,2,3,nres)
3596 ! write (iout,*) "gxcart(nres-1)", gloc_sc(2,nres-3,icg),
3597 ! & dtauangle(j,2,3,nres), gxcart(j,nres-1)
3598 ! if (itype(nres-2,1).ne.10) then
3599 if ((itype(nres-2,1).ne.10).and.&
3600 (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2))).and.(molnum(nres-2).ne.5)) then
3601 gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) &
3602 *dtauangle(j,3,3,nres)
3604 if ((itype(nres,1).ne.10).and.&
3605 (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))).and.(molnum(nres).ne.5)) then
3606 gxcart(j,nres-1)=gxcart(j,nres-1)-gloc_sc(3,nres-2,icg) &
3607 *dtauangle(j,3,1,nres+1)
3608 gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(3,nres-2,icg) &
3609 *dtauangle(j,3,2,nres+1)
3612 if ((itype(nres-2,1).ne.10).and.&
3613 (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2))).and.(molnum(nres-2).ne.5)) then
3614 gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(1,nres-3,icg)* &
3615 dtauangle(j,1,3,nres)
3617 if ((itype(nres,1).ne.10).and.(itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))).and.(molnum(nres).ne.5)) then
3618 gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(2,nres-2,icg)* &
3619 dtauangle(j,2,2,nres+1)
3620 ! write (iout,*) "gcart(nres-1)", gloc_sc(2,nres-2,icg),
3621 ! & dtauangle(j,2,2,nres+1), itype(nres-1,1),itype(nres,1)
3626 if ((nres.ge.3).and.(itype(nres,1).ne.10).and. &
3627 (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))).and.(molnum(nres).ne.5))then
3629 gxcart(j,nres)=gxcart(j,nres)+gloc_sc(3,nres-2,icg) &
3630 *dtauangle(j,3,3,nres+1)+gloc_sc(2,nres-2,icg) &
3631 *dtauangle(j,2,3,nres+1)
3634 ! write(iout,*) "final gcart",gcart(1,2)
3635 ! The side-chain vector derivatives
3636 ! print *,"gcart",gcart(:,:)
3638 end subroutine int_to_cart
3639 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3640 !-----------------------------------------------------------------------------
3642 !-----------------------------------------------------------------------------
3643 subroutine gen_dist_constr
3644 ! Generate CA distance constraints.
3645 ! implicit real*8 (a-h,o-z)
3646 ! include 'DIMENSIONS'
3647 ! include 'COMMON.IOUNITS'
3648 ! include 'COMMON.GEO'
3649 ! include 'COMMON.VAR'
3650 ! include 'COMMON.INTERACT'
3651 ! include 'COMMON.LOCAL'
3652 ! include 'COMMON.NAMES'
3653 ! include 'COMMON.CHAIN'
3654 ! include 'COMMON.FFIELD'
3655 ! include 'COMMON.SBRIDGE'
3656 ! include 'COMMON.HEADER'
3657 ! include 'COMMON.CONTROL'
3658 ! include 'COMMON.DBASE'
3659 ! include 'COMMON.THREAD'
3660 ! include 'COMMON.TIME1'
3661 ! integer :: itype_pdb !(maxres)
3662 ! common /pizda/ itype_pdb(nres)
3663 character(len=2) :: iden
3666 !d print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
3667 !d write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
3668 !d & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
3670 do i=nstart_sup,nstart_sup+nsup-1
3671 !d write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
3672 !d & ' seq_pdb', restyp(itype_pdb(i))
3673 do j=i+2,nstart_sup+nsup-1
3675 ihpb(nhpb)=i+nstart_seq-nstart_sup
3676 jhpb(nhpb)=j+nstart_seq-nstart_sup
3678 dhpb(nhpb)=dist(i,j)
3681 !d write (iout,'(a)') 'Distance constraints:'
3686 !d if (ii.gt.nres) then
3691 !d write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
3692 !d & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
3693 !d & dhpb(i),forcon(i)
3695 ! deallocate(itype_pdb)
3698 end subroutine gen_dist_constr
3700 !-----------------------------------------------------------------------------
3702 !-----------------------------------------------------------------------------
3703 subroutine cartprint
3705 use geometry_data, only: c
3706 use energy_data, only: itype
3707 ! implicit real*8 (a-h,o-z)
3708 ! include 'DIMENSIONS'
3709 ! include 'COMMON.CHAIN'
3710 ! include 'COMMON.INTERACT'
3711 ! include 'COMMON.NAMES'
3712 ! include 'COMMON.IOUNITS'
3717 write (iout,110) restyp(itype(i,1),1),i,c(1,i),c(2,i),&
3718 c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
3720 100 format (//' alpha-carbon coordinates ',&
3721 ' centroid coordinates'/ &
3722 ' ', 6X,'X',11X,'Y',11X,'Z',&
3723 10X,'X',11X,'Y',11X,'Z')
3724 110 format (a,'(',i3,')',6f12.5)
3726 end subroutine cartprint
3727 !-----------------------------------------------------------------------------
3728 !-----------------------------------------------------------------------------
3729 subroutine alloc_geo_arrays
3730 !EL Allocation of tables used by module energy
3732 integer :: i,j,nres2
3736 allocate(phibound(2,nres+2)) !(2,maxres)
3737 !----------------------
3739 ! common /chain/ in molread
3740 ! real(kind=8),dimension(:,:),allocatable :: c !(3,maxres2+2)
3741 ! real(kind=8),dimension(:,:),allocatable :: dc
3742 allocate(dc_old(3,0:nres2))
3743 ! if(.not.allocated(dc_norm2)) allocate(dc_norm2(3,0:nres2+2)) !(3,0:maxres2)
3744 if(.not.allocated(dc_norm2)) then
3745 allocate(dc_norm2(3,0:nres2+2)) !(3,0:maxres2)
3749 !el if(.not.allocated(dc_norm))
3750 !elwrite(iout,*) "jestem w alloc geo 1"
3751 if(.not.allocated(dc_norm)) then
3752 allocate(dc_norm(3,0:nres2+2)) !(3,0:maxres2)
3755 !elwrite(iout,*) "jestem w alloc geo 1"
3756 allocate(xloc(3,nres),xrot(3,nres))
3757 !elwrite(iout,*) "jestem w alloc geo 1"
3759 !elwrite(iout,*) "jestem w alloc geo 1"
3760 allocate(dc_work(6*nres)) !(MAXRES6) maxres6=6*maxres
3762 allocate(t(3,3,nres),r(3,3,nres))
3763 allocate(prod(3,3,nres),rt(3,3,nres)) !(3,3,maxres)
3764 ! common /refstruct/
3765 if(.not.allocated(cref)) allocate(cref(3,nres2+2,maxperm)) !(3,maxres2+2,maxperm)
3766 !elwrite(iout,*) "jestem w alloc geo 2"
3767 allocate(crefjlee(3,nres2+2)) !(3,maxres2+2)
3768 if(.not.allocated(chain_rep)) allocate(chain_rep(3,nres2+2,maxsym)) !(3,maxres2+2,maxsym)
3769 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3770 ! common /from_zscore/ in module.compare
3771 !----------------------
3773 ! Inverses of the actual virtual bond lengths
3774 ! common /invlen/ in io_conf: molread or readpdb
3775 ! real(kind=8),dimension(:),allocatable :: vbld_inv !(maxres2)
3776 !----------------------
3778 ! Store the geometric variables in the following COMMON block.
3779 ! common /var/ in readpdb or ...
3780 if(.not.allocated(theta)) allocate(theta(nres+2))
3781 if(.not.allocated(phi)) allocate(phi(nres+2))
3782 if(.not.allocated(alph)) allocate(alph(nres+2))
3783 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3784 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3785 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3786 if(.not.allocated(costtab)) allocate(costtab(nres))
3787 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3788 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3789 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3790 ! real(kind=8),dimension(:),allocatable :: vbld !(2*maxres) in io_conf: molread or readpdb
3791 allocate(omicron(2,nres+2)) !(2,maxres)
3792 allocate(tauangle(3,nres+2)) !(3,maxres)
3793 !elwrite(iout,*) "jestem w alloc geo 3"
3794 if(.not.allocated(xxtab)) allocate(xxtab(nres))
3795 if(.not.allocated(yytab)) allocate(yytab(nres))
3796 if(.not.allocated(zztab)) allocate(zztab(nres)) !(maxres)
3797 if(.not.allocated(xxref)) allocate(xxref(nres))
3798 if(.not.allocated(yyref)) allocate(yyref(nres))
3799 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3800 allocate(ialph(nres,2)) !(maxres,2)
3803 allocate(ivar(4*nres2)) !(4*maxres2)
3805 #if defined(WHAM_RUN) || defined(CLUSTER)
3806 allocate(vbld(2*nres))
3808 allocate(vbld_inv(2*nres))
3813 end subroutine alloc_geo_arrays
3814 !-----------------------------------------------------------------------------
3815 !-----------------------------------------------------------------------------
3816 subroutine returnbox
3817 integer :: allareout,i,j,k,nojumpval,chain_beg,mnum
3818 integer :: chain_end,ireturnval
3819 real*8 :: difference
3820 !C change suggested by Ana - end
3824 !C write(*,*) 'initial', i,j,c(j,i)
3826 !C change suggested by Ana - begin
3828 !C change suggested by Ana -end
3831 if ((itype(i,mnum).eq.ntyp1_molec(mnum))&
3832 .and.(itype(i+1,mnum).eq.ntyp1_molec(mnum))) then
3834 if (allareout.eq.1) then
3835 ireturnval=int(c(j,i)/boxxsize)
3836 if (c(j,i).le.0) ireturnval=ireturnval-1
3837 do k=chain_beg,chain_end
3838 c(j,k)=c(j,k)-ireturnval*boxxsize
3839 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
3842 if (chain_beg.eq.1) &
3843 dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
3844 !C Suggested by Ana -end
3849 if (int(c(j,i)/boxxsize).eq.0) allareout=0
3852 if (allareout.eq.1) then
3853 ireturnval=int(c(j,i)/boxxsize)
3854 if (c(j,i).le.0) ireturnval=ireturnval-1
3856 c(j,k)=c(j,k)-ireturnval*boxxsize
3857 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
3862 !C write(*,*) 'befor no jump', i,j,c(j,i)
3867 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
3868 .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then
3869 difference=abs(c(j,i-1)-c(j,i))
3870 !C print *,'diff', difference
3871 if (difference.gt.boxxsize/2.0) then
3872 if (c(j,i-1).gt.c(j,i)) then
3881 c(j,i)=c(j,i)+nojumpval*boxxsize
3882 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
3887 if (itype(i,mnum).eq.ntyp1_molec(mnum) .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then
3888 difference=abs(c(j,i-1)-c(j,i))
3889 if (difference.gt.boxxsize/2.0) then
3890 if (c(j,i-1).gt.c(j,i)) then
3899 c(j,i)=c(j,i)+nojumpval*boxxsize
3900 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
3904 !C write(*,*) 'after no jump', i,j,c(j,i)
3908 !C suggesed by Ana begins
3914 if ((itype(i,mnum).eq.ntyp1_molec(mnum))&
3915 .and.(itype(i+1,mnum).eq.ntyp1_molec(mnum))) then
3917 if (allareout.eq.1) then
3918 ireturnval=int(c(j,i)/boxysize)
3919 if (c(j,i).le.0) ireturnval=ireturnval-1
3920 do k=chain_beg,chain_end
3921 c(j,k)=c(j,k)-ireturnval*boxysize
3922 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
3925 if (chain_beg.eq.1) &
3926 dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
3927 !C Suggested by Ana -end
3932 if (int(c(j,i)/boxysize).eq.0) allareout=0
3935 if (allareout.eq.1) then
3936 ireturnval=int(c(j,i)/boxysize)
3937 if (c(j,i).le.0) ireturnval=ireturnval-1
3939 c(j,k)=c(j,k)-ireturnval*boxysize
3940 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
3946 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
3947 .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then
3948 difference=abs(c(j,i-1)-c(j,i))
3949 if (difference.gt.boxysize/2.0) then
3950 if (c(j,i-1).gt.c(j,i)) then
3959 c(j,i)=c(j,i)+nojumpval*boxysize
3960 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
3965 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
3966 .and. itype(i-1,mnum).eq.ntyp1) then
3967 difference=abs(c(j,i-1)-c(j,i))
3968 if (difference.gt.boxysize/2.0) then
3969 if (c(j,i-1).gt.c(j,i)) then
3978 c(j,i)=c(j,i)+nojumpval*boxysize
3979 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
3982 !C Suggested by Ana -begins
3984 !C Suggested by Ana -ends
3989 if ((itype(i,mnum).eq.ntyp1_molec(mnum))&
3990 .and.(itype(i+1,mnum).eq.ntyp1_molec(mnum))) then
3992 if (allareout.eq.1) then
3993 ireturnval=int(c(j,i)/boxysize)
3994 if (c(j,i).le.0) ireturnval=ireturnval-1
3995 do k=chain_beg,chain_end
3996 c(j,k)=c(j,k)-ireturnval*boxzsize
3997 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
4000 if (chain_beg.eq.1) dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
4001 !C Suggested by Ana -end
4006 if (int(c(j,i)/boxzsize).eq.0) allareout=0
4009 if (allareout.eq.1) then
4010 ireturnval=int(c(j,i)/boxzsize)
4011 if (c(j,i).le.0) ireturnval=ireturnval-1
4013 c(j,k)=c(j,k)-ireturnval*boxzsize
4014 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
4020 if (itype(i,mnum).eq.ntyp1_molec(mnum) .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then
4021 difference=abs(c(j,i-1)-c(j,i))
4022 if (difference.gt.(boxzsize/2.0)) then
4023 if (c(j,i-1).gt.c(j,i)) then
4032 c(j,i)=c(j,i)+nojumpval*boxzsize
4033 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
4038 if (itype(i,mnum).eq.ntyp1_molec(mnum) &
4039 .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then
4040 difference=abs(c(j,i-1)-c(j,i))
4041 if (difference.gt.boxzsize/2.0) then
4042 if (c(j,i-1).gt.c(j,i)) then
4051 c(j,i)=c(j,i)+nojumpval*boxzsize
4052 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
4055 if (molnum(i).eq.5) then
4056 c(1,i)=dmod(c(1,i),boxxsize)
4057 if (c(1,i).lt.0) c(1,i)=c(1,i)+boxxsize
4058 c(2,i)=dmod(c(2,i),boxysize)
4059 if (c(2,i).lt.0) c(2,i)=c(2,i)+boxysize
4060 c(3,i)=dmod(c(3,i),boxzsize)
4061 if (c(3,i).lt.0) c(3,i)=c(3,i)+boxzsize
4062 c(1,i+nres)=dmod(c(1,i+nres),boxxsize)
4063 c(2,i+nres)=dmod(c(2,i+nres),boxysize)
4064 c(3,i+nres)=dmod(c(3,i+nres),boxzsize)
4068 end subroutine returnbox
4069 !-------------------------------------------------------------------------------------------------------