1 subroutine gen_rand_conf(nstart,nend,*)
2 C Generate random conformation or chain cut and regrowth.
3 implicit real*8 (a-h,o-z)
8 include 'COMMON.INTERACT'
9 include 'COMMON.IOUNITS'
12 include 'COMMON.CONTROL'
13 logical overlap,back,fail,ldir
15 write (iout,*) ' CG Processor',me,' maxgen=',maxgen
18 write (iout,*) 'Gen_Rand_conf: nstart=',nstart,' nend',nend,
25 phi(4)=gen_phi(4,iabs(itype(2)),iabs(itype(3)))
26 c write(iout,*)'phi(4)=',rad2deg*phi(4)
27 if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2)),pi,phi(4))
28 c write(iout,*)'theta(3)=',rad2deg*theta(3)
32 do while (fail.and.nsi.le.maxsi)
33 call gen_side(it1,theta(3),alph(2),omeg(2),fail)
36 if (nsi.gt.maxsi) return1
51 do while (i.le.nend .and. niter.lt.maxgen)
54 write (iout,'(/80(1h*)/2a/80(1h*))')
55 & 'Generation procedure went down to ',
56 & 'chain beginning. Cannot continue...'
57 write (*,'(/80(1h*)/2a/80(1h*))')
58 & 'Generation procedure went down to ',
59 & 'chain beginning. Cannot continue...'
66 c print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
67 c & ' nit=',nit,' niter=',niter,' maxgen=',maxgen
68 phi(i+1)=gen_phi(i+1,it1,it)
70 phi(i)=gen_phi(i+1,it2,it1)
71 c print *,'phi(',i,')=',phi(i)
72 theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
73 if (it2.ne.10 .and. it2.ne.ntyp1) then
76 do while (fail.and.nsi.le.maxsi)
77 call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail)
80 if (nsi.gt.maxsi) return1
82 call locate_next_res(i-1)
84 theta(i)=gen_theta(it1,phi(i),phi(i+1))
85 if (it1.ne.10 .and. it1.ne.ntyp1) then
88 do while (fail.and.nsi.le.maxsi)
89 call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
92 if (nsi.gt.maxsi) return1
94 call locate_next_res(i)
95 if (overlap(i-1,ldir)) then
96 if (nit.lt.maxnit) then
106 & 'Cannot generate non-overlaping conformation. Increase MAXNIT.'
108 & 'Cannot generate non-overlaping conformation. Increase MAXNIT.'
119 if (niter.ge.maxgen) then
120 write (iout,'(a,2i5)')
121 & 'Too many trials in conformation generation',niter,maxgen
123 & 'Too many trials in conformation generation',niter,maxgen
128 c(j,nres+nres)=c(j,nres)
139 do while (i.ge.nend .and. niter.lt.maxgen)
140 if (i.gt.nstart) then
142 write (iout,'(/80(1h*)/2a/80(1h*))')
143 & 'Generation procedure went down to ',
144 & 'chain beginning. Cannot continue...'
145 write (*,'(/80(1h*)/2a/80(1h*))')
146 & 'Generation procedure went down to ',
147 & 'chain beginning. Cannot continue...'
154 c print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
155 c & ' nit=',nit,' niter=',niter,' maxgen=',maxgen
156 phi(i+2)=gen_phi(i+2,it1,it2)
158 phi(i+3)=gen_phi(i+3,it2,it3)
159 phi(i+4)=gen_phi(i+3,it2,it3)
160 c print *,'phi(',i+3,')=',phi(i+3)
161 theta(i+3)=gen_theta(it2,phi(i+3),phi(i+4))
162 if (it2.ne.10 .and. it2.ne.ntyp1) then
165 do while (fail.and.nsi.le.maxsi)
166 call gen_side(it2,theta(i+3),alph(i+2),omeg(i+2),fail)
169 if (nsi.gt.maxsi) return1
171 call locate_prev_res(i+1)
173 theta(i+2)=gen_theta(it1,phi(i+2),phi(i+3))
174 write (iout,*) "i"," theta",theta(i+2)," phi",phi(i+2),phi(i+3)
175 if (it1.ne.10 .and. it1.ne.ntyp1) then
178 do while (fail.and.nsi.le.maxsi)
179 call gen_side(it1,theta(i+2),alph(i+1),omeg(i+1),fail)
182 write (iout,*) "i",i," alpha",alph(i+1)," omeg",omeg(i+1),
184 if (nsi.gt.maxsi) return1
186 call locate_prev_res(i)
187 write (iout,*) "After locate_prev_res i=",i
188 write (iout,'(3f10.5,5x,3f10.5)') (c(l,i),l=1,3),
189 & (c(l,i+nres),l=1,3)
190 write (iout,'(3f10.5,5x,3f10.5)') (c(l,i+1),l=1,3),
191 & (c(l,i+1+nres),l=1,3)
193 if (overlap(i+1,ldir)) then
194 write (iout,*) "OVERLAP",nit
196 if (nit.lt.maxnit) then
206 & 'Cannot generate non-overlaping conformation. Increase MAXNIT.'
208 & 'Cannot generate non-overlaping conformation. Increase MAXNIT.'
213 write (iout,*) "================>> NO OVERLAP",i
220 write (iout,*) "iter",niter
222 if (niter.ge.maxgen) then
223 write (iout,'(a,2i5)')
224 &'Too many trials in backward conformation generation',niter,maxgen
226 &'Too many trials in backward conformation generation',niter,maxgen
232 write (iout,*) "HERE!!!!"
236 c-------------------------------------------------------------------------
237 logical function overlap(i,ldir)
238 implicit real*8 (a-h,o-z)
240 include 'COMMON.IOUNITS'
241 include 'COMMON.CHAIN'
242 include 'COMMON.INTERACT'
243 include 'COMMON.FFIELD'
246 logical lprn /.false./
249 if (iti.gt.ntyp) return
250 C Check for SC-SC overlaps.
251 cd print *,'nnt=',nnt,' nct=',nct
252 c write (iout,*) "overlap i",i," ldir",ldir
258 if (j.lt.i-1 .or. ipot.ne.4) then
259 rcomp=sigmaii(iti,itj)
264 if (dist(nres+i,nres+j).lt.redfac*rcomp) then
266 if (lprn) write (iout,*) 'overlap, SC-SC: i=',i,' j=',j,
267 & ' dist=',dist(nres+i,nres+j),' rcomp=',
272 C Check for overlaps between the added peptide group and the preceding
276 c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
280 if (dist(nres+j,maxres2+1).lt.4.0D0*redfac) then
281 if (lprn) write (iout,*) 'overlap, p-Sc: i=',i,' j=',j,
282 & ' dist=',dist(nres+j,maxres2+1)
287 C Check for overlaps between the added side chain and the preceding peptide
291 c(k,maxres2+1)=0.5D0*(c(k,j)+c(k,j+1))
293 if (dist(nres+i,maxres2+1).lt.4.0D0*redfac) then
294 if (lprn) write (iout,*) 'overlap, SC-p: i=',i,' j=',j,
295 & ' dist=',dist(nres+i,maxres2+1)
300 C Check for p-p overlaps
302 c(j,maxres2+2)=0.5D0*(c(j,i)+c(j,i+1))
307 c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1))
309 if(iteli.ne.0.and.itelj.ne.0)then
310 if (dist(maxres2+1,maxres2+2).lt.rpp(iteli,itelj)*redfac) then
311 if (lprn) write (iout,*) 'overlap, p-p: i=',i,' j=',j,
312 & ' dist=',dist(maxres2+1,maxres2+2)
321 if (lprn) write (iout,*) "start",i+1," end",nres_start+nsup
322 do j=i+1,nres_start+nsup
324 if (j.gt.i+1 .or. ipot.ne.4) then
325 rcomp=sigmaii(iti,itj)
330 c write (iout,*) "SCSC j",j," dist",dist(nres+i,nres+j),
331 c & " radfac",redfac*rcomp
332 if (dist(nres+i,nres+j).lt.redfac*rcomp) then
334 if (lprn) write (iout,*) 'overlap, SC-SC: i=',i,' j=',j,
335 & ' dist=',dist(nres+i,nres+j),' rcomp=',
340 C Check for overlaps between the added peptide group and the succeeding
344 c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
346 do j=i+2,nstart_res+nsup
348 cd write (iout,*) "pSC j",j," dist",dist(nres+j,maxres2+1)
349 if (dist(nres+j,maxres2+1).lt.4.0D0*redfac) then
350 if (lprn) write (iout,*) 'overlap, p-Sc: i=',i,' j=',j,
351 & ' dist=',dist(nres+j,maxres2+1)
356 C Check for overlaps between the added side chain and the succeeding peptide
358 do j=i+2,nstart_seq+nsup
360 c(k,maxres2+1)=0.5D0*(c(k,j)+c(k,j+1))
362 cd write (iout,*) "SCp j",j," dist",dist(nres+i,maxres2+1)
363 if (dist(nres+i,maxres2+1).lt.4.0D0*redfac) then
364 if (lprn) write (iout,*) 'overlap, SC-p: i=',i,' j=',j,
365 & ' dist=',dist(nres+i,maxres2+1)
370 C Check for p-p overlaps
372 c(j,maxres2+2)=0.5D0*(c(j,i)+c(j,i+1))
374 do j=i+2,nres_start+nsup
377 c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1))
379 if(iteli.ne.0.and.itelj.ne.0)then
380 cd write (iout,*) "pp j",j," dist",dist(maxres2+1,maxres2+2)
381 if (dist(maxres2+1,maxres2+2).lt.rpp(iteli,itelj)*redfac) then
382 if (lprn) write (iout,*) 'overlap, p-p: i=',i,' j=',j,
383 & ' dist=',dist(maxres2+1,maxres2+2)
394 c--------------------------------------------------------------------------
395 double precision function gen_phi(i,it1,it2)
396 implicit real*8 (a-h,o-z)
398 include "COMMON.TORCNSTR"
400 include 'COMMON.BOUNDS'
401 if (raw_psipred .or. ndih_constr.eq.0) then
402 gen_phi=ran_number(-pi,pi)
404 C 8/13/98 Generate phi using pre-defined boundaries
405 gen_phi=ran_number(phibound(1,i),phibound(2,i))
409 c---------------------------------------------------------------------------
410 double precision function gen_theta(it,gama,gama1)
411 implicit real*8 (a-h,o-z)
413 include 'COMMON.LOCAL'
415 double precision y(2),z(2)
416 double precision theta_max,theta_min
417 c print *,'gen_theta: it=',it
420 if (dabs(gama).gt.dwapi) then
427 if (dabs(gama1).gt.dwapi) then
434 thet_pred_mean=a0thet(it)
436 thet_pred_mean=thet_pred_mean+athet(k,it,1,1)*y(k)
437 & +bthet(k,it,1,1)*z(k)
441 sig=sig*thet_pred_mean+polthet(j,it)
443 sig=0.5D0/(sig*sig+sigc0(it))
445 &0.5D0*((gthet(2,it)-thet_pred_mean)/gthet(3,it))**2)
446 c print '(i5,5(1pe14.4))',it,(gthet(j,it),j=1,3)
447 c print '(5(1pe14.4))',thet_pred_mean,theta0(it),sig,sig0(it),ak
448 theta_temp=binorm(thet_pred_mean,theta0(it),sig,sig0(it),ak)
449 if (theta_temp.lt.theta_min) theta_temp=theta_min
450 if (theta_temp.gt.theta_max) theta_temp=theta_max
452 c print '(a)','Exiting GENTHETA.'
455 c-------------------------------------------------------------------------
456 subroutine gen_side(it,the,al,om,fail)
457 implicit real*8 (a-h,o-z)
460 include 'COMMON.LOCAL'
461 include 'COMMON.SETUP'
462 include 'COMMON.IOUNITS'
463 double precision MaxBoxLen /10.0D0/
464 double precision Ap_inv(3,3),a(3,3),z(3,maxlob),W1(maxlob),
465 & sumW(0:maxlob),y(2),cm(2),eig(2),box(2,2),work(100),detAp(maxlob)
466 double precision eig_limit /1.0D-8/
467 double precision Big /10.0D0/
468 double precision vec(3,3)
469 logical lprint,fail,lcheck
473 if (the.eq.0.0D0 .or. the.eq.pi) then
475 write (*,'(a,i4,a,i3,a,1pe14.5)')
476 & 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the
478 cd write (iout,'(a,i3,a,1pe14.5)')
479 cd & 'Error in GenSide: it=',it,' theta=',the
488 print '(a,i4,a)','CG Processor:',me,' Enter Gen_Side.'
489 write (iout,'(a,i4,a)') 'Processor:',me,' Enter Gen_Side.'
491 print *,'it=',it,' nlobit=',nlobit,' the=',the,' tant=',tant
492 write (iout,*) 'it=',it,' nlobit=',nlobit,' the=',the,
496 zz1=tant-censc(1,i,it)
499 a(k,l)=gaussc(k,l,i,it)
502 detApi=a(2,2)*a(3,3)-a(2,3)**2
503 Ap_inv(2,2)=a(3,3)/detApi
504 Ap_inv(2,3)=-a(2,3)/detApi
505 Ap_inv(3,2)=Ap_inv(2,3)
506 Ap_inv(3,3)=a(2,2)/detApi
508 write (*,'(/a,i2/)') 'Cluster #',i
509 write (*,'(3(1pe14.5),5x,1pe14.5)')
510 & ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
511 write (iout,'(/a,i2/)') 'Cluster #',i
512 write (iout,'(3(1pe14.5),5x,1pe14.5)')
513 & ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
518 W1i=W1i+a(k,1)*a(l,1)*Ap_inv(k,l)
522 W1(i)=dexp(bsc(i,it)-0.5D0*W1i*zz1*zz1)
523 c if (lprint) write(*,'(a,3(1pe15.5)/)')
524 c & 'detAp, W1, anormi',detApi,W1i,anormi
528 zk=zk+zz1*Ap_inv(k,l)*a(l,1)
532 detAp(i)=dsqrt(detApi)
536 print *,'W1:',(w1(i),i=1,nlobit)
537 print *,'detAp:',(detAp(i),i=1,nlobit)
540 print '(i2,3f10.5)',i,(rad2deg*z(j,i),j=2,3)
542 write (iout,*) 'W1:',(w1(i),i=1,nlobit)
543 write (iout,*) 'detAp:',(detAp(i),i=1,nlobit)
546 write (iout,'(i2,3f10.5)') i,(rad2deg*z(j,i),j=2,3)
550 C Writing the distribution just to check the procedure
556 fac=fac+W1(i)/detAp(i)
558 fac=1.0D0/(2.0D0*fac*pi)
559 cd print *,it,'fac=',fac
568 a(j-1,k-1)=gaussc(j,k,i,it)
580 wykl=wykl+a(j,k)*(y(j)-z(j+1,i))*(y(k)-z(k+1,i))
583 wart=wart+W1(i)*dexp(-0.5D0*wykl)
590 c print *,'y',y(1),y(2),' fac=',fac
592 write (20,'(2f10.3,1pd15.5)') y(1)*rad2deg,y(2)*rad2deg,wart
597 c print *,'it=',it,' sum=',sum*dV,' sum1=',sum1*dV
601 C Calculate the CM of the system
608 sumW(i)=sumW(i-1)+W1(i)
613 cm(1)=cm(1)+z(2,j)*W1(j)
614 cm(2)=cm(2)+W1(j)*(z(3,1)+pinorm(z(3,j)-z(3,1)))
616 cm(1)=cm(1)/sumW(nlobit)
617 cm(2)=cm(2)/sumW(nlobit)
618 if (cm(1).gt.Big .or. cm(1).lt.-Big .or.
619 & cm(2).gt.Big .or. cm(2).lt.-Big) then
620 cd write (iout,'(a)')
621 cd & 'Unexpected error in GenSide - CM coordinates too large.'
622 cd write (iout,'(i5,2(1pe14.5))') it,cm(1),cm(2)
624 cd & 'Unexpected error in GenSide - CM coordinates too large.'
625 cd write (*,'(i5,2(1pe14.5))') it,cm(1),cm(2)
629 cd print *,'CM:',cm(1),cm(2)
631 C Find the largest search distance from CM
637 a(j-1,k-1)=gaussc(j,k,i,it)
641 call f02faf('N','U',2,a,3,eig,work,100,ifail)
643 call djacob(2,3,10000,1.0d-10,a,vec,eig)
647 print *,'*************** CG Processor',me
648 print *,'CM:',cm(1),cm(2)
649 write (iout,*) '*************** CG Processor',me
650 write (iout,*) 'CM:',cm(1),cm(2)
651 print '(A,8f10.5)','Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
652 write (iout,'(A,8f10.5)')
653 & 'Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
656 if (eig(1).lt.eig_limit) then
658 & 'From Mult_Norm: Eigenvalues of A are too small.'
660 & 'From Mult_Norm: Eigenvalues of A are too small.'
667 radius=radius+pinorm(z(j+1,i)-cm(j))**2
669 radius=dsqrt(radius)+3.0D0/dsqrt(eig(1))
670 if (radius.gt.radmax) radmax=radius
672 if (radmax.gt.pi) radmax=pi
674 C Determine the boundaries of the search rectangle.
677 print '(a,4(1pe14.4))','W1: ',(W1(i),i=1,nlob(it) )
678 print '(a,4(1pe14.4))','radmax: ',radmax
680 box(1,1)=dmax1(cm(1)-radmax,0.0D0)
681 box(2,1)=dmin1(cm(1)+radmax,pi)
682 box(1,2)=cm(2)-radmax
683 box(2,2)=cm(2)+radmax
686 print *,'CG Processor',me,' Array BOX:'
690 print '(4(1pe14.4))',((box(k,j),k=1,2),j=1,2)
691 print '(a,4(1pe14.4))','sumW: ',(sumW(i),i=0,nlob(it) )
693 write (iout,*)'CG Processor',me,' Array BOX:'
695 write (iout,*)'Array BOX:'
697 write(iout,'(4(1pe14.4))') ((box(k,j),k=1,2),j=1,2)
698 write(iout,'(a,4(1pe14.4))')'sumW: ',(sumW(i),i=0,nlob(it) )
700 if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
702 write (iout,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
703 write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
705 c write (iout,'(a)') 'Bad sampling box.'
710 which_lobe=ran_number(0.0D0,sumW(nlobit))
711 c print '(a,1pe14.4)','which_lobe=',which_lobe
713 if (sumW(i-1).le.which_lobe .and. sumW(i).ge.which_lobe) goto 1
716 c print *,'ilob=',ilob,' nlob=',nlob(it)
720 a(i-1,j-1)=gaussc(i,j,ilob,it)
723 cd print '(a,i4,a)','CG Processor',me,' Calling MultNorm1.'
724 call mult_norm1(3,2,a,cm,box,y,fail)
728 cd print *,'al=',al,' om=',om
732 c---------------------------------------------------------------------------
733 double precision function ran_number(x1,x2)
734 C Calculate a random real number from the range (x1,x2).
735 implicit real*8 (a-h,o-z)
737 double precision x1,x2,fctor
738 data fctor /2147483647.0D0/
741 include 'COMMON.SETUP'
742 ran_number=x1+(x2-x1)*prng_next(me)
745 ran_number=x1+(x2-x1)*ix/fctor
749 c--------------------------------------------------------------------------
750 integer function iran_num(n1,n2)
751 C Calculate a random integer number from the range (n1,n2).
752 implicit real*8 (a-h,o-z)
755 real fctor /2147483647.0/
758 include 'COMMON.SETUP'
759 ix=n1+(n2-n1+1)*prng_next(me)
765 ix=n1+(n2-n1+1)*(ix/fctor)
771 c--------------------------------------------------------------------------
772 double precision function binorm(x1,x2,sigma1,sigma2,ak)
773 implicit real*8 (a-h,o-z)
774 c print '(a)','Enter BINORM.'
775 alowb=dmin1(x1-3.0D0*sigma1,x2-3.0D0*sigma2)
776 aupb=dmax1(x1+3.0D0*sigma1,x2+3.0D0*sigma2)
777 seg=sigma1/(sigma1+ak*sigma2)
778 alen=ran_number(0.0D0,1.0D0)
779 if (alen.lt.seg) then
780 binorm=anorm_distr(x1,sigma1,alowb,aupb)
782 binorm=anorm_distr(x2,sigma2,alowb,aupb)
784 c print '(a)','Exiting BINORM.'
787 c-----------------------------------------------------------------------
788 c double precision function anorm_distr(x,sigma,alowb,aupb)
789 c implicit real*8 (a-h,o-z)
790 c print '(a)','Enter ANORM_DISTR.'
791 c 10 y=ran_number(alowb,aupb)
792 c expon=dexp(-0.5D0*((y-x)/sigma)**2)
793 c ran=ran_number(0.0D0,1.0D0)
794 c if (expon.lt.ran) goto 10
796 c print '(a)','Exiting ANORM_DISTR.'
799 c-----------------------------------------------------------------------
800 double precision function anorm_distr(x,sigma,alowb,aupb)
801 implicit real*8 (a-h,o-z)
802 c to make a normally distributed deviate with zero mean and unit variance
805 real fac,gset,rsq,v1,v2,ran1
809 1 v1=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
810 v2=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
812 if(rsq.ge.1.d0.or.rsq.eq.0.0d0) goto 1
813 fac=sqrt(-2.0d0*log(rsq)/rsq)
821 anorm_distr=x+gaussdev*sigma
824 c------------------------------------------------------------------------
825 subroutine mult_norm(lda,n,a,x,fail)
827 C Generate the vector X whose elements obey the multiple-normal distribution
828 C from exp(-0.5*X'AX). LDA is the leading dimension of the moment matrix A,
829 C n is the dimension of the problem. FAIL is set at .TRUE., if the smallest
830 C eigenvalue of the matrix A is close to 0.
832 implicit double precision (a-h,o-z)
833 double precision a(lda,n),x(n),eig(100),vec(3,3),work(100)
834 double precision eig_limit /1.0D-8/
837 c print '(a)','Enter MULT_NORM.'
839 C Find the smallest eigenvalue of the matrix A.
842 c print '(8f10.5)',(a(i,j),j=1,n)
845 call f02faf('V','U',2,a,lda,eig,work,100,ifail)
847 call djacob(2,lda,10000,1.0d-10,a,vec,eig)
849 c print '(8f10.5)',(eig(i),i=1,n)
852 c print '(8f10.5)',(a(i,j),j=1,n)
854 if (eig(1).lt.eig_limit) then
855 print *,'From Mult_Norm: Eigenvalues of A are too small.'
860 C Generate points following the normal distributions along the principal
861 C axes of the moment matrix. Store in WORK.
864 sigma=1.0D0/dsqrt(eig(i))
866 work(i)=anorm_distr(0.0D0,sigma,-alim,alim)
869 C Transform the vector of normal variables back to the original basis.
880 c------------------------------------------------------------------------
881 subroutine mult_norm1(lda,n,a,z,box,x,fail)
883 C Generate the vector X whose elements obey the multi-gaussian multi-dimensional
884 C distribution from sum_{i=1}^m W(i)exp[-0.5*X'(i)A(i)X(i)]. LDA is the
885 C leading dimension of the moment matrix A, n is the dimension of the
886 C distribution, nlob is the number of lobes. FAIL is set at .TRUE., if the
887 C smallest eigenvalue of the matrix A is close to 0.
889 implicit real*8 (a-h,o-z)
894 double precision a(lda,n),z(n),x(n),box(n,n)
895 double precision etmp
896 include 'COMMON.IOUNITS'
898 include 'COMMON.SETUP'
902 C Generate points following the normal distributions along the principal
903 C axes of the moment matrix. Store in WORK.
905 cd print *,'CG Processor',me,' entered MultNorm1.'
906 cd print '(2(1pe14.4),3x,1pe14.4)',((a(i,j),j=1,2),z(i),i=1,2)
908 cd print *,i,box(1,i),box(2,i)
912 if (istep.gt.10000) then
913 c write (iout,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
915 c write (*,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
917 c write (iout,*) 'box',box
918 c write (iout,*) 'a',a
919 c write (iout,*) 'z',z
924 x(i)=ran_number(box(1,i),box(2,i))
929 ww=ww+0.5D0*a(i,i)*xi*xi
931 ww=ww+a(i,j)*xi*pinorm(x(j)-z(j))
934 dec=ran_number(0.0D0,1.0D0)
935 c print *,(x(i),i=1,n),ww,dexp(-ww),dec
936 crc if (dec.gt.dexp(-ww)) goto 10
942 if (dec.gt.etmp) goto 10
943 cd print *,'CG Processor',me,' exitting MultNorm1.'
947 crc--------------------------------------
948 subroutine overlap_sc(scfail)
949 c Internal and cartesian coordinates must be consistent as input,
950 c and will be up-to-date on return.
951 c At the end of this procedure, scfail is true if there are
952 c overlapping residues left, or false otherwise (success)
953 implicit real*8 (a-h,o-z)
955 include 'COMMON.CHAIN'
956 include 'COMMON.INTERACT'
957 include 'COMMON.FFIELD'
959 include 'COMMON.SBRIDGE'
960 include 'COMMON.IOUNITS'
961 logical had_overlaps,fail,scfail
962 integer ioverlap(maxres),ioverlap_last
965 call overlap_sc_list(ioverlap,ioverlap_last)
966 if (ioverlap_last.gt.0) then
967 write (iout,*) '#OVERLAPing residues ',ioverlap_last
968 write (iout,'(20i4)') (ioverlap(k),k=1,ioverlap_last)
974 if (ioverlap_last.eq.0) exit
976 do ires=1,ioverlap_last
982 do while (fail.and.nsi.le.maxsi)
983 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
990 call chainbuild_extconf
991 call overlap_sc_list(ioverlap,ioverlap_last)
992 c write (iout,*) 'Overlaping residues ',ioverlap_last,
993 c & (ioverlap(j),j=1,ioverlap_last)
996 if (k.le.1000.and.ioverlap_last.eq.0) then
998 if (had_overlaps) then
999 write (iout,*) '#OVERLAPing all corrected after ',k,
1000 & ' random generation'
1004 write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
1005 write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last)
1011 write (iout,'(a30,i5,a12,i4)')
1012 & '#OVERLAP FAIL in gen_side after',maxsi,
1018 subroutine overlap_sc_list(ioverlap,ioverlap_last)
1019 implicit real*8 (a-h,o-z)
1020 include 'DIMENSIONS'
1021 include 'COMMON.GEO'
1022 include 'COMMON.LOCAL'
1023 include 'COMMON.IOUNITS'
1024 include 'COMMON.CHAIN'
1025 include 'COMMON.INTERACT'
1026 include 'COMMON.FFIELD'
1027 include 'COMMON.VAR'
1028 include 'COMMON.CALC'
1030 integer ioverlap(maxres),ioverlap_last
1034 C Check for SC-SC overlaps and mark residues
1035 c print *,'>>overlap_sc nnt=',nnt,' nct=',nct
1037 do i=iatsc_s,iatsc_e
1038 itypi=iabs(itype(i))
1039 itypi1=iabs(itype(i+1))
1043 dxi=dc_norm(1,nres+i)
1044 dyi=dc_norm(2,nres+i)
1045 dzi=dc_norm(3,nres+i)
1046 dsci_inv=dsc_inv(itypi)
1048 do iint=1,nint_gr(i)
1049 do j=istart(i,iint),iend(i,iint)
1051 itypj=iabs(itype(j))
1052 dscj_inv=dsc_inv(itypj)
1053 sig0ij=sigma(itypi,itypj)
1054 chi1=chi(itypi,itypj)
1055 chi2=chi(itypj,itypi)
1062 alf12=0.5D0*(alf1+alf2)
1064 rcomp=sigmaii(itypi,itypj)
1066 rcomp=sigma(itypi,itypj)
1068 c print '(2(a3,2i3),a3,2f10.5)',
1069 c & ' i=',i,iti,' j=',j,itj,' d=',dist(nres+i,nres+j)
1074 dxj=dc_norm(1,nres+j)
1075 dyj=dc_norm(2,nres+j)
1076 dzj=dc_norm(3,nres+j)
1077 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1081 sig=sig0ij*dsqrt(sigsq)
1082 rij_shift=1.0D0/rij-sig+sig0ij
1084 ct if ( 1.0/rij .lt. redfac*rcomp .or.
1085 ct & rij_shift.le.0.0D0 ) then
1086 if ( rij_shift.le.0.0D0 ) then
1087 cd write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)')
1088 cd & 'overlap SC-SC: i=',i,' j=',j,
1089 cd & ' dist=',dist(nres+i,nres+j),' rcomp=',
1090 cd & rcomp,1.0/rij,rij_shift
1091 ioverlap_last=ioverlap_last+1
1092 ioverlap(ioverlap_last)=i
1093 do k=1,ioverlap_last-1
1094 if (ioverlap(k).eq.i) ioverlap_last=ioverlap_last-1
1096 ioverlap_last=ioverlap_last+1
1097 ioverlap(ioverlap_last)=j
1098 do k=1,ioverlap_last-1
1099 if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1