1 subroutine gen_rand_conf(nstart,*)
2 C Generate random conformation or chain cut and regrowth.
8 include 'COMMON.INTERACT'
9 include 'COMMON.IOUNITS'
12 include 'COMMON.CONTROL'
13 logical overlap,back,fail
15 integer i,j,k,it,it1,it2,nit,niter,nsi,maxsi,maxnit
16 double precision gen_theta,gen_phi,dist,ran_number
17 cd print *,' CG Processor',me,' maxgen=',maxgen
19 cd write (iout,*) 'Gen_Rand_conf: nstart=',nstart
22 phi(4)=gen_phi(4,iabs(itype(2)),iabs(itype(3)))
23 c write(iout,*)'phi(4)=',rad2deg*phi(4)
24 if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2)),pi,phi(4))
25 c write(iout,*)'theta(3)=',rad2deg*theta(3)
29 do while (fail.and.nsi.le.maxsi)
30 call gen_side(it1,theta(3),alph(2),omeg(2),fail)
33 if (nsi.gt.maxsi) return1
48 do while (i.le.nres .and. niter.lt.maxgen)
51 write (iout,'(/80(1h*)/2a/80(1h*))')
52 & 'Generation procedure went down to ',
53 & 'chain beginning. Cannot continue...'
54 write (*,'(/80(1h*)/2a/80(1h*))')
55 & 'Generation procedure went down to ',
56 & 'chain beginning. Cannot continue...'
63 if (it.eq.ntyp1 .and. it1.eq.ntyp1) then
64 vbld(i)=ran_number(3.8d0,10.0d0)
65 vbld_inv(i)=1.0d0/vbld(i)
67 c write (iout,*) 'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,
68 c & ' it2=',it2,' nit=',nit,' niter=',niter,' maxgen=',maxgen
69 phi(i+1)=gen_phi(i+1,it1,it)
71 phi(i)=gen_phi(i+1,it2,it1)
72 c write(iout,*) 'phi(',i,')=',phi(i)," type",it1,it2,it
73 theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
74 if (it2.ne.10 .and. it2.ne.ntyp1) then
77 do while (fail.and.nsi.le.maxsi)
78 call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail)
81 if (nsi.gt.maxsi) return1
83 call locate_next_res(i-1)
85 if (it1.ne.ntyp1) then
86 theta(i)=gen_theta(it1,phi(i),phi(i+1))
88 theta(i)=ran_number(1.326d0,2.548d0)
90 c write (iout,*) "i",i," theta",theta(i)
91 if (it1.ne.10 .and. it1.ne.ntyp1) then
94 do while (fail.and.nsi.le.maxsi)
95 call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
98 c write (iout,*) "After forward SC generation:",nsi,maxsi
99 if (nsi.gt.maxsi) return1
101 call locate_next_res(i)
102 if (overlap(i-1)) then
103 c write (iout,*) "overlap",i-1
104 if (nit.lt.maxnit) then
114 & 'Cannot generate non-overlaping conformation. Increase MAXNIT.'
116 & 'Cannot generate non-overlaping conformation. Increase MAXNIT.'
121 c write (iout,*) "No overlap",i-1
128 if (niter.ge.maxgen) then
129 write (iout,'(a,2i5)')
130 & 'Too many trials in conformation generation',niter,maxgen
132 & 'Too many trials in conformation generation',niter,maxgen
137 c(j,nres+nres)=c(j,nres)
141 c-------------------------------------------------------------------------
142 logical function overlap(i)
145 include 'COMMON.CHAIN'
146 include 'COMMON.INTERACT'
147 include 'COMMON.FFIELD'
148 double precision redfac /0.5D0/
149 integer i,j,k,iti,itj,iteli,itelj
150 double precision rcomp
151 double precision dist
154 if (iti.gt.ntyp) return
155 C Check for SC-SC overlaps.
156 cd print *,'nnt=',nnt,' nct=',nct
159 if (itj.eq.ntyp1) cycle
160 if (j.lt.i-1 .or. ipot.ne.4) then
161 rcomp=sigmaii(iti,itj)
166 if (dist(nres+i,nres+j).lt.redfac*rcomp) then
168 c print *,'overlap, SC-SC: i=',i,' j=',j,
169 c & ' dist=',dist(nres+i,nres+j),' rcomp=',
174 C Check for overlaps between the added peptide group and the preceding
178 c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
182 cd print *,'overlap, p-Sc: i=',i,' j=',j,
183 cd & ' dist=',dist(nres+j,maxres2+1)
184 if (dist(nres+j,maxres2+1).lt.4.0D0*redfac) then
189 C Check for overlaps between the added side chain and the preceding peptide
193 c(k,maxres2+1)=0.5D0*(c(k,j)+c(k,j+1))
195 cd print *,'overlap, SC-p: i=',i,' j=',j,
196 cd & ' dist=',dist(nres+i,maxres2+1)
197 if (dist(nres+i,maxres2+1).lt.4.0D0*redfac) then
202 C Check for p-p overlaps
204 c(j,maxres2+2)=0.5D0*(c(j,i)+c(j,i+1))
209 c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1))
211 cd print *,'overlap, p-p: i=',i,' j=',j,
212 cd & ' dist=',dist(maxres2+1,maxres2+2)
213 if(iteli.ne.0.and.itelj.ne.0)then
214 if (dist(maxres2+1,maxres2+2).lt.rpp(iteli,itelj)*redfac) then
222 c--------------------------------------------------------------------------
223 double precision function gen_phi(i,it1,it2)
224 implicit real*8 (a-h,o-z)
226 include 'COMMON.IOUNITS'
227 include "COMMON.TORCNSTR"
229 include 'COMMON.BOUNDS'
230 include 'COMMON.INTERACT'
231 double precision sumprob(3)
232 double precision pinorm
234 if (ndih_constr.eq.0) then
235 gen_phi=ran_number(-pi,pi)
236 else if (raw_psipred) then
237 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
238 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
240 sumprob(1)=vpsipred(2,ii)
241 sumprob(2)=sumprob(1)+vpsipred(3,ii)
242 sumprob(3)=sumprob(2)+vpsipred(1,ii)
243 aux=ran_number(0.0d0,sumprob(3))
245 write(iout,*)"gen_phi: residue i",i," ii",ii," vpsipred",
246 & vpsipred(2,ii),vpsipred(3,ii),vpsipred(1,ii)," sumprob",
247 & sumprob(1),sumprob(2),sumprob(3)
248 write (iout,*) "aux",aux
250 if (aux.le.sumprob(1)) then
252 write (iout,*) "hel:",
253 & (phibound(1,i)-sdihed(1,ndih_constr))*rad2deg,
254 & (phibound(1,i)+sdihed(1,ndih_constr))*rad2deg
256 gen_phi=ran_number(phibound(1,i)-sdihed(1,ndih_constr)
257 & ,phibound(1,i)+sdihed(1,ndih_constr))
258 else if (aux.le.sumprob(2)) then
260 write (iout,*) "ext:",
261 & (phibound(2,i)-sdihed(2,ndih_constr))*rad2deg,
262 & (phibound(2,i)+sdihed(2,ndih_constr))*rad2deg
264 gen_phi=pinorm(ran_number(phibound(2,i)-sdihed(2,ndih_constr)
265 & ,phibound(2,i)+sdihed(2,ndih_constr)))
268 write (iout,*) "coil:",-180.0,180.0
270 gen_phi=ran_number(-pi,pi)
273 gen_phi=ran_number(-pi,pi)
276 C 8/13/98 Generate phi using pre-defined boundaries
277 gen_phi=ran_number(phibound(1,i),phibound(2,i))
281 c---------------------------------------------------------------------------
282 double precision function gen_theta(it,gama,gama1)
283 implicit real*8 (a-h,o-z)
285 include 'COMMON.LOCAL'
287 double precision y(2),z(2)
288 double precision theta_max,theta_min
289 c print *,'gen_theta: it=',it
292 if (dabs(gama).gt.dwapi) then
299 if (dabs(gama1).gt.dwapi) then
306 thet_pred_mean=a0thet(it)
308 thet_pred_mean=thet_pred_mean+athet(k,it,1,1)*y(k)
309 & +bthet(k,it,1,1)*z(k)
313 sig=sig*thet_pred_mean+polthet(j,it)
315 sig=0.5D0/(sig*sig+sigc0(it))
317 &0.5D0*((gthet(2,it)-thet_pred_mean)/gthet(3,it))**2)
318 c print '(i5,5(1pe14.4))',it,(gthet(j,it),j=1,3)
319 c print '(5(1pe14.4))',thet_pred_mean,theta0(it),sig,sig0(it),ak
320 theta_temp=binorm(thet_pred_mean,theta0(it),sig,sig0(it),ak)
321 if (theta_temp.lt.theta_min) theta_temp=theta_min
322 if (theta_temp.gt.theta_max) theta_temp=theta_max
324 c print '(a)','Exiting GENTHETA.'
327 c-------------------------------------------------------------------------
328 subroutine gen_side(it,the,al,om,fail)
329 implicit real*8 (a-h,o-z)
332 include 'COMMON.LOCAL'
333 include 'COMMON.SETUP'
334 include 'COMMON.IOUNITS'
335 double precision MaxBoxLen /10.0D0/
336 double precision Ap_inv(3,3),a(3,3),z(3,maxlob),W1(maxlob),
337 & sumW(0:maxlob),y(2),cm(2),eig(2),box(2,2),work(100),detAp(maxlob)
338 double precision eig_limit /1.0D-8/
339 double precision Big /10.0D0/
340 double precision vec(3,3)
341 logical lprint,fail,lcheck,lprn /.false./
345 if (the.eq.0.0D0 .or. the.eq.pi) then
347 write (iout,'(a,i4,a,i3,a,1pe14.5)')
348 & 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the
350 cd write (iout,'(a,i3,a,1pe14.5)')
351 cd & 'Error in GenSide: it=',it,' theta=',the
360 print '(a,i4,a)','CG Processor:',me,' Enter Gen_Side.'
361 write (iout,'(a,i4,a)') 'Processor:',me,' Enter Gen_Side.'
363 print *,'it=',it,' nlobit=',nlobit,' the=',the,' tant=',tant
364 write (iout,*) 'it=',it,' nlobit=',nlobit,' the=',the,
368 zz1=tant-censc(1,i,it)
371 a(k,l)=gaussc(k,l,i,it)
374 detApi=a(2,2)*a(3,3)-a(2,3)**2
375 Ap_inv(2,2)=a(3,3)/detApi
376 Ap_inv(2,3)=-a(2,3)/detApi
377 Ap_inv(3,2)=Ap_inv(2,3)
378 Ap_inv(3,3)=a(2,2)/detApi
380 write (*,'(/a,i2/)') 'Cluster #',i
381 write (*,'(3(1pe14.5),5x,1pe14.5)')
382 & ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
383 write (iout,'(/a,i2/)') 'Cluster #',i
384 write (iout,'(3(1pe14.5),5x,1pe14.5)')
385 & ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
390 W1i=W1i+a(k,1)*a(l,1)*Ap_inv(k,l)
394 W1(i)=dexp(bsc(i,it)-0.5D0*W1i*zz1*zz1)
395 c if (lprint) write(*,'(a,3(1pe15.5)/)')
396 c & 'detAp, W1, anormi',detApi,W1i,anormi
400 zk=zk+zz1*Ap_inv(k,l)*a(l,1)
404 detAp(i)=dsqrt(detApi)
408 print *,'W1:',(w1(i),i=1,nlobit)
409 print *,'detAp:',(detAp(i),i=1,nlobit)
412 print '(i2,3f10.5)',i,(rad2deg*z(j,i),j=2,3)
414 write (iout,*) 'W1:',(w1(i),i=1,nlobit)
415 write (iout,*) 'detAp:',(detAp(i),i=1,nlobit)
418 write (iout,'(i2,3f10.5)') i,(rad2deg*z(j,i),j=2,3)
422 C Writing the distribution just to check the procedure
428 fac=fac+W1(i)/detAp(i)
430 fac=1.0D0/(2.0D0*fac*pi)
431 cd print *,it,'fac=',fac
440 a(j-1,k-1)=gaussc(j,k,i,it)
452 wykl=wykl+a(j,k)*(y(j)-z(j+1,i))*(y(k)-z(k+1,i))
455 wart=wart+W1(i)*dexp(-0.5D0*wykl)
462 c print *,'y',y(1),y(2),' fac=',fac
464 write (20,'(2f10.3,1pd15.5)') y(1)*rad2deg,y(2)*rad2deg,wart
469 c print *,'it=',it,' sum=',sum*dV,' sum1=',sum1*dV
473 C Calculate the CM of the system
480 sumW(i)=sumW(i-1)+W1(i)
485 cm(1)=cm(1)+z(2,j)*W1(j)
486 cm(2)=cm(2)+W1(j)*(z(3,1)+pinorm(z(3,j)-z(3,1)))
488 cm(1)=cm(1)/sumW(nlobit)
489 cm(2)=cm(2)/sumW(nlobit)
490 if (cm(1).gt.Big .or. cm(1).lt.-Big .or.
491 & cm(2).gt.Big .or. cm(2).lt.-Big) then
492 cd write (iout,'(a)')
493 cd & 'Unexpected error in GenSide - CM coordinates too large.'
494 cd write (iout,'(i5,2(1pe14.5))') it,cm(1),cm(2)
496 cd & 'Unexpected error in GenSide - CM coordinates too large.'
497 cd write (*,'(i5,2(1pe14.5))') it,cm(1),cm(2)
501 cd print *,'CM:',cm(1),cm(2)
503 C Find the largest search distance from CM
509 a(j-1,k-1)=gaussc(j,k,i,it)
513 call f02faf('N','U',2,a,3,eig,work,100,ifail)
515 call djacob(2,3,10000,1.0d-10,a,vec,eig)
519 print *,'*************** CG Processor',me
520 print *,'CM:',cm(1),cm(2)
521 write (iout,*) '*************** CG Processor',me
522 write (iout,*) 'CM:',cm(1),cm(2)
523 print '(A,8f10.5)','Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
524 write (iout,'(A,8f10.5)')
525 & 'Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
528 if (eig(1).lt.eig_limit) then
530 & 'From Mult_Norm: Eigenvalues of A are too small.'
532 & 'From Mult_Norm: Eigenvalues of A are too small.'
539 radius=radius+pinorm(z(j+1,i)-cm(j))**2
541 radius=dsqrt(radius)+3.0D0/dsqrt(eig(1))
542 if (radius.gt.radmax) radmax=radius
544 if (radmax.gt.pi) radmax=pi
546 C Determine the boundaries of the search rectangle.
549 print '(a,4(1pe14.4))','W1: ',(W1(i),i=1,nlob(it) )
550 print '(a,4(1pe14.4))','radmax: ',radmax
552 box(1,1)=dmax1(cm(1)-radmax,0.0D0)
553 box(2,1)=dmin1(cm(1)+radmax,pi)
554 box(1,2)=cm(2)-radmax
555 box(2,2)=cm(2)+radmax
558 print *,'CG Processor',me,' Array BOX:'
562 print '(4(1pe14.4))',((box(k,j),k=1,2),j=1,2)
563 print '(a,4(1pe14.4))','sumW: ',(sumW(i),i=0,nlob(it) )
565 write (iout,*)'CG Processor',me,' Array BOX:'
567 write (iout,*)'Array BOX:'
569 write(iout,'(4(1pe14.4))') ((box(k,j),k=1,2),j=1,2)
570 write(iout,'(a,4(1pe14.4))')'sumW: ',(sumW(i),i=0,nlob(it) )
572 if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
575 write (iout,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
576 write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
579 if (lprn) write (iout,'(a)') 'Bad sampling box.'
584 which_lobe=ran_number(0.0D0,sumW(nlobit))
585 c print '(a,1pe14.4)','which_lobe=',which_lobe
587 if (sumW(i-1).le.which_lobe .and. sumW(i).ge.which_lobe) goto 1
590 c print *,'ilob=',ilob,' nlob=',nlob(it)
594 a(i-1,j-1)=gaussc(i,j,ilob,it)
597 cd print '(a,i4,a)','CG Processor',me,' Calling MultNorm1.'
598 call mult_norm1(3,2,a,cm,box,y,fail)
602 cd print *,'al=',al,' om=',om
606 c---------------------------------------------------------------------------
607 double precision function ran_number(x1,x2)
608 C Calculate a random real number from the range (x1,x2).
609 implicit real*8 (a-h,o-z)
611 double precision x1,x2,fctor
612 data fctor /2147483647.0D0/
615 include 'COMMON.SETUP'
616 ran_number=x1+(x2-x1)*prng_next(me)
619 ran_number=x1+(x2-x1)*ix/fctor
623 c--------------------------------------------------------------------------
624 integer function iran_num(n1,n2)
625 C Calculate a random integer number from the range (n1,n2).
626 implicit real*8 (a-h,o-z)
629 real fctor /2147483647.0/
632 include 'COMMON.SETUP'
633 ix=n1+(n2-n1+1)*prng_next(me)
639 ix=n1+(n2-n1+1)*(ix/fctor)
645 c--------------------------------------------------------------------------
646 double precision function binorm(x1,x2,sigma1,sigma2,ak)
647 implicit real*8 (a-h,o-z)
648 c print '(a)','Enter BINORM.'
649 alowb=dmin1(x1-3.0D0*sigma1,x2-3.0D0*sigma2)
650 aupb=dmax1(x1+3.0D0*sigma1,x2+3.0D0*sigma2)
651 seg=sigma1/(sigma1+ak*sigma2)
652 alen=ran_number(0.0D0,1.0D0)
653 if (alen.lt.seg) then
654 binorm=anorm_distr(x1,sigma1,alowb,aupb)
656 binorm=anorm_distr(x2,sigma2,alowb,aupb)
658 c print '(a)','Exiting BINORM.'
661 c-----------------------------------------------------------------------
662 c double precision function anorm_distr(x,sigma,alowb,aupb)
663 c implicit real*8 (a-h,o-z)
664 c print '(a)','Enter ANORM_DISTR.'
665 c 10 y=ran_number(alowb,aupb)
666 c expon=dexp(-0.5D0*((y-x)/sigma)**2)
667 c ran=ran_number(0.0D0,1.0D0)
668 c if (expon.lt.ran) goto 10
670 c print '(a)','Exiting ANORM_DISTR.'
673 c-----------------------------------------------------------------------
674 double precision function anorm_distr(x,sigma,alowb,aupb)
675 implicit real*8 (a-h,o-z)
676 c to make a normally distributed deviate with zero mean and unit variance
679 real fac,gset,rsq,v1,v2,ran1
683 1 v1=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
684 v2=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
686 if(rsq.ge.1.d0.or.rsq.eq.0.0d0) goto 1
687 fac=sqrt(-2.0d0*log(rsq)/rsq)
695 anorm_distr=x+gaussdev*sigma
698 c------------------------------------------------------------------------
699 subroutine mult_norm(lda,n,a,x,fail)
701 C Generate the vector X whose elements obey the multiple-normal distribution
702 C from exp(-0.5*X'AX). LDA is the leading dimension of the moment matrix A,
703 C n is the dimension of the problem. FAIL is set at .TRUE., if the smallest
704 C eigenvalue of the matrix A is close to 0.
706 implicit double precision (a-h,o-z)
707 double precision a(lda,n),x(n),eig(100),vec(3,3),work(100)
708 double precision eig_limit /1.0D-8/
711 c print '(a)','Enter MULT_NORM.'
713 C Find the smallest eigenvalue of the matrix A.
716 c print '(8f10.5)',(a(i,j),j=1,n)
719 call f02faf('V','U',2,a,lda,eig,work,100,ifail)
721 call djacob(2,lda,10000,1.0d-10,a,vec,eig)
723 c print '(8f10.5)',(eig(i),i=1,n)
726 c print '(8f10.5)',(a(i,j),j=1,n)
728 if (eig(1).lt.eig_limit) then
729 print *,'From Mult_Norm: Eigenvalues of A are too small.'
734 C Generate points following the normal distributions along the principal
735 C axes of the moment matrix. Store in WORK.
738 sigma=1.0D0/dsqrt(eig(i))
740 work(i)=anorm_distr(0.0D0,sigma,-alim,alim)
743 C Transform the vector of normal variables back to the original basis.
754 c------------------------------------------------------------------------
755 subroutine mult_norm1(lda,n,a,z,box,x,fail)
757 C Generate the vector X whose elements obey the multi-gaussian multi-dimensional
758 C distribution from sum_{i=1}^m W(i)exp[-0.5*X'(i)A(i)X(i)]. LDA is the
759 C leading dimension of the moment matrix A, n is the dimension of the
760 C distribution, nlob is the number of lobes. FAIL is set at .TRUE., if the
761 C smallest eigenvalue of the matrix A is close to 0.
763 implicit real*8 (a-h,o-z)
768 double precision a(lda,n),z(n),x(n),box(n,n)
769 double precision etmp
770 include 'COMMON.IOUNITS'
772 include 'COMMON.SETUP'
776 C Generate points following the normal distributions along the principal
777 C axes of the moment matrix. Store in WORK.
779 cd print *,'CG Processor',me,' entered MultNorm1.'
780 cd print '(2(1pe14.4),3x,1pe14.4)',((a(i,j),j=1,2),z(i),i=1,2)
782 cd print *,i,box(1,i),box(2,i)
786 if (istep.gt.10000) then
787 c write (iout,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
789 c write (*,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
791 c write (iout,*) 'box',box
792 c write (iout,*) 'a',a
793 c write (iout,*) 'z',z
798 x(i)=ran_number(box(1,i),box(2,i))
803 ww=ww+0.5D0*a(i,i)*xi*xi
805 ww=ww+a(i,j)*xi*pinorm(x(j)-z(j))
808 dec=ran_number(0.0D0,1.0D0)
809 c print *,(x(i),i=1,n),ww,dexp(-ww),dec
810 crc if (dec.gt.dexp(-ww)) goto 10
816 if (dec.gt.etmp) goto 10
817 cd print *,'CG Processor',me,' exitting MultNorm1.'
821 crc--------------------------------------
822 subroutine overlap_sc(scfail)
823 c Internal and cartesian coordinates must be consistent as input,
824 c and will be up-to-date on return.
825 c At the end of this procedure, scfail is true if there are
826 c overlapping residues left, or false otherwise (success)
827 implicit real*8 (a-h,o-z)
829 include 'COMMON.CHAIN'
830 include 'COMMON.INTERACT'
831 include 'COMMON.FFIELD'
833 include 'COMMON.SBRIDGE'
834 include 'COMMON.IOUNITS'
835 logical had_overlaps,fail,scfail
836 integer ioverlap(maxres),ioverlap_last
837 integer maxit_corr /5000/
840 call overlap_sc_list(ioverlap,ioverlap_last,.false.)
841 if (ioverlap_last.gt.0) then
842 write (iout,*) '#OVERLAPing residues ',ioverlap_last
843 write (iout,'(15i6)') (ioverlap(k),k=1,ioverlap_last)
849 if (ioverlap_last.eq.0) exit
851 do ires=1,ioverlap_last
854 if (iti.ne.10 .and. iti.lt.ntyp1) then
857 do while (fail.and.nsi.le.maxsi)
858 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
859 call sc_coord_rebuild(i)
865 c write (iout,*) "before chaincuild overlap_sc_list: dc0",dc(:,0)
866 c call chainbuild_extconf
867 c write (iout,*) "after chaincuild overlap_sc_list: dc0",dc(:,0)
868 call overlap_sc_list(ioverlap,ioverlap_last,.false.)
869 write (iout,*)'#Overlaping residues @iter',k,":",ioverlap_last
870 write (iout,*)'Residue list:',(ioverlap(j),j=1,ioverlap_last)
873 if (k.le.maxit_corr.and.ioverlap_last.eq.0) then
875 if (had_overlaps) then
876 write (iout,*) '#OVERLAPing all corrected after ',k,
877 & ' random generation'
881 write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
882 write (iout,'(15i6)') (ioverlap(j),j=1,ioverlap_last)
888 write (iout,'(a30,i5,a12,i6)')
889 & '#OVERLAP FAIL in gen_side after',maxsi,
895 subroutine overlap_sc_list(ioverlap,ioverlap_last,lprn)
900 include "COMMON.SETUP"
902 integer ioverlap_last_tab(0:max_fg_procs-1),
903 & ioverlap_all(maxres*max_fg_procs),displs(0:max_fg_procs-1)
906 include 'COMMON.LOCAL'
907 include 'COMMON.IOUNITS'
908 include 'COMMON.CHAIN'
909 include 'COMMON.INTERACT'
910 include 'COMMON.FFIELD'
912 include 'COMMON.CALC'
913 integer ii,itypi,itypj,itypi1,ind,ikont
915 integer ioverlap(maxres),ioverlap_last
916 double precision redfac /0.5D0/
917 double precision rrij,rij_shift,sig0ij,xi,yi,zi,rcomp,sig
918 double precision dist
921 if (nfgtasks.gt.1) then
923 & call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
924 call MPI_Bcast(c(1,1),6*nres,MPI_DOUBLE_PRECISION,king,FG_COMM,
926 call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,king,
928 call MPI_Bcast(dc_norm(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,
929 & king,FG_COMM,IERROR)
932 c write (iout,*) "overlap_sc_list"
933 c write(iout,*) "iatsc_s",iatsc_s," iatsc_e",iatsc_e
934 c write(iout,*) "nnt",nnt," nct",nct
936 C Check for SC-SC overlaps and mark residues
937 c print *,'>>overlap_sc nnt=',nnt,' nct=',nct
940 write (iout,*) "FG proecssor",fg_rank," g_listscsc_start",
941 & g_listscsc_start," g_listscsc_end",g_listscsc_end
942 write (*,*) "FG proecssor",fg_rank," g_listscsc_start",
943 & g_listscsc_start," g_listscsc_end",g_listscsc_end
945 c do i=iatsc_s,iatsc_e
946 do ikont=g_listscsc_start,g_listscsc_end
947 c write (*,*) "FG processor",fg_rank," loop begins ioverlap_last",
949 i=newcontlisti(ikont)
950 j=newcontlistj(ikont)
952 c write (*,*) "FG processor",fg_rank," loop begins ioverlap_last",
953 c & ioverlap_last,"ikont i j",ikont,i,j
955 itypi1=iabs(itype(i+1))
956 if (itypi.eq.ntyp1) cycle
960 dxi=dc_norm(1,nres+i)
961 dyi=dc_norm(2,nres+i)
962 dzi=dc_norm(3,nres+i)
963 dsci_inv=dsc_inv(itypi)
965 c do iint=1,nint_gr(i)
966 c do j=istart(i,iint),iend(i,iint)
970 if (itypj.eq.ntyp1) cycle
971 c write (iout,*) "i,j",i,j," itypi,itypj",itypi,itypj
972 dscj_inv=dsc_inv(itypj)
973 sig0ij=sigma(itypi,itypj)
974 chi1=chi(itypi,itypj)
975 chi2=chi(itypj,itypi)
982 alf12=0.5D0*(alf1+alf2)
984 rcomp=sigmaii(itypi,itypj)
986 rcomp=sigma(itypi,itypj)
988 c write (iout,'(2(a3,2i5),a3,2f10.5)'),
989 c & ' i=',i,itypi,' j=',j,itypj,' d=',dist(nres+i,nres+j)
994 dxj=dc_norm(1,nres+j)
995 dyj=dc_norm(2,nres+j)
996 dzj=dc_norm(3,nres+j)
997 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1000 c write (iout,*) "dxj",dxj," dyj",dyj," dzj",dzj
1001 c write (iout,*) "erij",erij
1002 c write (iout,*) "om1",om1," om2",om2," om12",om12,
1003 c & " faceps1",faceps1," eps1",eps1
1004 c write (iout,*) "sigsq",sigsq
1006 sig=sig0ij*dsqrt(sigsq)
1007 rij_shift=1.0D0/rij-sig+sig0ij
1008 c write (iout,*) "rij",1.0d0/rij," rij_shift",rij_shift,
1009 c & " sig",sig," sig0ij",sig0ij
1010 c if ( rij_shift.le.0.0D0 ) then
1011 if ( rij_shift/sig0ij.le.0.1D0 ) then
1012 c write (iout,*) "overlap",i,j
1014 write (iout,'(a,i5,a,i5,a,f10.5,a,3f10.5)')
1015 & 'overlap SC-SC: i=',i,' j=',j,
1016 & ' dist=',dist(nres+i,nres+j),' rcomp=',
1017 & rcomp,1.0/rij,rij_shift
1018 write (*,'(a,i2,a,i5,a,i5,a,f10.5,a,3f10.5)')
1019 & 'FG processor',fg_rank,' overlap SC-SC: i=',i,' j=',j,
1020 & ' dist=',dist(nres+i,nres+j),' rcomp=',
1021 & rcomp,1.0/rij,rij_shift
1023 ioverlap_last=ioverlap_last+1
1024 ioverlap(ioverlap_last)=i
1025 do k=1,ioverlap_last-1
1026 if (ioverlap(k).eq.i) ioverlap_last=ioverlap_last-1
1028 ioverlap_last=ioverlap_last+1
1029 ioverlap(ioverlap_last)=j
1030 do k=1,ioverlap_last-1
1031 if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1
1033 c write(*,*) "FG processor",fg_rank,i,j," ioverlap_last",
1034 c & ioverlap_last," ioverlap",(ioverlap(k),k=1,ioverlap_last)
1041 write (iout,*) "FG Processor",fg_rank," ioverlap_last",
1042 & ioverlap_last," ioverlap",(ioverlap(i),i=1,ioverlap_last)
1043 write (*,*) "FG Processor",fg_rank," ioverlap_last",ioverlap_last,
1044 & " ioverlap",(ioverlap(i),i=1,ioverlap_last)
1047 if (nfgtasks.eq.1) return
1049 write (iout,*) "Before MPI_Gather"
1052 call MPI_Gather(ioverlap_last,1,MPI_INTEGER,ioverlap_last_tab,
1053 & 1,MPI_INTEGER,king,FG_COMM,IERROR)
1055 write (iout,*) "After MPI_Gather"
1060 & write (iout,*) "FG Processor",fg_rank,"ioverlap_last_tab",
1061 & (ioverlap_last_tab(i),i=0,nfgtasks-1)
1066 displs(i)=displs(i-1)+ioverlap_last_tab(i-1)
1068 call MPI_Gatherv(ioverlap,ioverlap_last,MPI_INTEGER,
1069 & ioverlap_all,ioverlap_last_tab,displs,MPI_INTEGER,king,
1072 write (iout,*) "After Gatherv"
1075 if (fg_rank.gt.0) return
1078 ioverlap_last=ioverlap_last+ioverlap_last_tab(i)
1081 write (iout,*) "ioverlap_last",ioverlap_last," ioverlap_last",
1082 & (ioverlap_all(i),i=1,ioverlap_last)
1086 do i=1,ioverlap_last
1087 ioverlap(ii+1)=ioverlap_all(i)
1089 if (ioverlap(ii+1).eq.ioverlap(j)) goto 11
1096 write (iout,*) "After summing: ioverlap_last",ioverlap_last,
1097 & " ioverlap",(ioverlap(i),i=1,ioverlap_last)