9be8b232d8a5ff085e8fd7bd227e32b16cf5e550
[unres.git] / source / unres / src_CSA_DiL / gen_rand_conf.F
1       subroutine gen_rand_conf(nstart,*)
2 C Generate random conformation or chain cut and regrowth.
3       implicit real*8 (a-h,o-z)
4       include 'DIMENSIONS'
5       include 'COMMON.CHAIN'
6       include 'COMMON.LOCAL'
7       include 'COMMON.VAR'
8       include 'COMMON.INTERACT'
9       include 'COMMON.IOUNITS'
10       include 'COMMON.MCM'
11       include 'COMMON.GEO'
12       include 'COMMON.CONTROL'
13       logical overlap,back,fail
14 cd    print *,' CG Processor',me,' maxgen=',maxgen
15       maxsi=100
16 cd    write (iout,*) 'Gen_Rand_conf: nstart=',nstart
17       if (nstart.lt.5) then
18         it1=itype(2)
19         phi(4)=gen_phi(4,itype(2),itype(3))
20 c       write(iout,*)'phi(4)=',rad2deg*phi(4)
21         if (nstart.lt.3) theta(3)=gen_theta(itype(2),pi,phi(4))
22 c       write(iout,*)'theta(3)=',rad2deg*theta(3) 
23         if (it1.ne.10) then
24           nsi=0
25           fail=.true.
26           do while (fail.and.nsi.le.maxsi)
27             call gen_side(it1,theta(3),alph(2),omeg(2),fail)
28             nsi=nsi+1
29           enddo
30           if (nsi.gt.maxsi) return1
31         endif ! it1.ne.10
32         call orig_frame
33         i=4
34         nstart=4
35       else
36         i=nstart
37         nstart=max0(i,4)
38       endif
39
40       maxnit=0
41
42       nit=0
43       niter=0
44       back=.false.
45       do while (i.le.nres .and. niter.lt.maxgen)
46         if (i.lt.nstart) then
47           if(iprint.gt.1) then
48           write (iout,'(/80(1h*)/2a/80(1h*))') 
49      &          'Generation procedure went down to ',
50      &          'chain beginning. Cannot continue...'
51           write (*,'(/80(1h*)/2a/80(1h*))') 
52      &          'Generation procedure went down to ',
53      &          'chain beginning. Cannot continue...'
54           endif
55           return1
56         endif
57         it1=itype(i-1)
58         it2=itype(i-2)
59         it=itype(i)
60 c       print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
61 c     &    ' nit=',nit,' niter=',niter,' maxgen=',maxgen
62         phi(i+1)=gen_phi(i+1,it1,it)
63         if (back) then
64           phi(i)=gen_phi(i+1,it2,it1)
65           print *,'phi(',i,')=',phi(i)
66           theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
67           if (it2.ne.10) then
68             nsi=0
69             fail=.true.
70             do while (fail.and.nsi.le.maxsi)
71               call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail)
72               nsi=nsi+1
73             enddo
74             if (nsi.gt.maxsi) return1
75           endif
76           call locate_next_res(i-1)
77         endif
78         theta(i)=gen_theta(it1,phi(i),phi(i+1))
79         if (it1.ne.10) then 
80         nsi=0
81         fail=.true.
82         do while (fail.and.nsi.le.maxsi)
83           call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
84           nsi=nsi+1
85         enddo
86         if (nsi.gt.maxsi) return1
87         endif
88         call locate_next_res(i)
89         if (overlap(i-1)) then
90           if (nit.lt.maxnit) then
91             back=.true.
92             nit=nit+1
93           else
94             nit=0
95             if (i.gt.3) then
96               back=.true.
97               i=i-1
98             else
99               write (iout,'(a)') 
100      &  'Cannot generate non-overlaping conformation. Increase MAXNIT.'
101               write (*,'(a)') 
102      &  'Cannot generate non-overlaping conformation. Increase MAXNIT.'
103               return1
104             endif
105           endif
106         else
107           back=.false.
108           nit=0 
109           i=i+1
110         endif
111         niter=niter+1
112       enddo
113       if (niter.ge.maxgen) then
114         write (iout,'(a,2i5)') 
115      & 'Too many trials in conformation generation',niter,maxgen
116         write (*,'(a,2i5)') 
117      & 'Too many trials in conformation generation',niter,maxgen
118         return1
119       endif
120       do j=1,3
121         c(j,nres+1)=c(j,1)
122         c(j,nres+nres)=c(j,nres)
123       enddo
124       return
125       end
126 c-------------------------------------------------------------------------
127       logical function overlap(i)
128       implicit real*8 (a-h,o-z)
129       include 'DIMENSIONS'
130       include 'COMMON.CHAIN'
131       include 'COMMON.INTERACT'
132       include 'COMMON.FFIELD'
133       data redfac /0.5D0/
134       overlap=.false.
135       iti=itype(i)
136       if (iti.gt.ntyp) return
137 C Check for SC-SC overlaps.
138 cd    print *,'nnt=',nnt,' nct=',nct
139       do j=nnt,i-1
140         itj=itype(j)
141         if (j.lt.i-1 .or. ipot.ne.4) then
142           rcomp=sigmaii(iti,itj)
143         else 
144           rcomp=sigma(iti,itj)
145         endif
146 cd      print *,'j=',j
147         if (dist(nres+i,nres+j).lt.redfac*rcomp) then
148           overlap=.true.
149 c        print *,'overlap, SC-SC: i=',i,' j=',j,
150 c     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
151 c     &     rcomp
152           return
153         endif
154       enddo
155 C Check for overlaps between the added peptide group and the preceding
156 C SCs.
157       iteli=itel(i)
158       do j=1,3
159         c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
160       enddo
161       do j=nnt,i-2
162         itj=itype(j)
163 cd      print *,'overlap, p-Sc: i=',i,' j=',j,
164 cd   &         ' dist=',dist(nres+j,maxres2+1)
165         if (dist(nres+j,maxres2+1).lt.4.0D0*redfac) then
166           overlap=.true.
167           return
168         endif
169       enddo
170 C Check for overlaps between the added side chain and the preceding peptide
171 C groups.
172       do j=1,nnt-2
173         do k=1,3
174           c(k,maxres2+1)=0.5D0*(c(k,j)+c(k,j+1))
175         enddo
176 cd      print *,'overlap, SC-p: i=',i,' j=',j,
177 cd   &         ' dist=',dist(nres+i,maxres2+1)
178         if (dist(nres+i,maxres2+1).lt.4.0D0*redfac) then
179           overlap=.true.
180           return
181         endif
182       enddo
183 C Check for p-p overlaps
184       do j=1,3
185         c(j,maxres2+2)=0.5D0*(c(j,i)+c(j,i+1))
186       enddo
187       do j=nnt,i-2
188         itelj=itel(j)
189         do k=1,3
190           c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1))
191         enddo
192 cd      print *,'overlap, p-p: i=',i,' j=',j,
193 cd   &         ' dist=',dist(maxres2+1,maxres2+2)
194         if(iteli.ne.0.and.itelj.ne.0)then
195         if (dist(maxres2+1,maxres2+2).lt.rpp(iteli,itelj)*redfac) then
196           overlap=.true.
197           return
198         endif
199         endif
200       enddo
201       return
202       end
203 c--------------------------------------------------------------------------
204       double precision function gen_phi(i,it1,it2)
205       implicit real*8 (a-h,o-z)
206       include 'DIMENSIONS'
207       include 'COMMON.GEO'
208       include 'COMMON.BOUNDS'
209 c      gen_phi=ran_number(-pi,pi) 
210 C 8/13/98 Generate phi using pre-defined boundaries
211       gen_phi=ran_number(phibound(1,i),phibound(2,i)) 
212       return
213       end
214 c---------------------------------------------------------------------------
215       double precision function gen_theta(it,gama,gama1)
216       implicit real*8 (a-h,o-z)
217       include 'DIMENSIONS'
218       include 'COMMON.LOCAL'
219       include 'COMMON.GEO'
220       double precision y(2),z(2)
221       double precision theta_max,theta_min
222 c     print *,'gen_theta: it=',it
223       theta_min=0.05D0*pi
224       theta_max=0.95D0*pi
225       if (dabs(gama).gt.dwapi) then
226         y(1)=dcos(gama)
227         y(2)=dsin(gama)
228       else
229         y(1)=0.0D0
230         y(2)=0.0D0
231       endif
232       if (dabs(gama1).gt.dwapi) then
233         z(1)=dcos(gama1)
234         z(2)=dsin(gama1)
235       else
236         z(1)=0.0D0
237         z(2)=0.0D0
238       endif  
239       thet_pred_mean=a0thet(it)
240       do k=1,2
241         thet_pred_mean=thet_pred_mean+athet(k,it)*y(k)+bthet(k,it)*z(k)
242       enddo
243       sig=polthet(3,it)
244       do j=2,0,-1
245         sig=sig*thet_pred_mean+polthet(j,it)
246       enddo
247       sig=0.5D0/(sig*sig+sigc0(it))
248       ak=dexp(gthet(1,it)-
249      &0.5D0*((gthet(2,it)-thet_pred_mean)/gthet(3,it))**2)
250 c     print '(i5,5(1pe14.4))',it,(gthet(j,it),j=1,3)
251 c     print '(5(1pe14.4))',thet_pred_mean,theta0(it),sig,sig0(it),ak
252       theta_temp=binorm(thet_pred_mean,theta0(it),sig,sig0(it),ak) 
253       if (theta_temp.lt.theta_min) theta_temp=theta_min
254       if (theta_temp.gt.theta_max) theta_temp=theta_max
255       gen_theta=theta_temp
256 c     print '(a)','Exiting GENTHETA.'
257       return
258       end
259 c-------------------------------------------------------------------------
260       subroutine gen_side(it,the,al,om,fail)
261       implicit real*8 (a-h,o-z)
262       include 'DIMENSIONS'
263       include 'COMMON.GEO'
264       include 'COMMON.LOCAL'
265       include 'COMMON.SETUP'
266       include 'COMMON.IOUNITS'
267       double precision MaxBoxLen /10.0D0/
268       double precision Ap_inv(3,3),a(3,3),z(3,maxlob),W1(maxlob),
269      & sumW(0:maxlob),y(2),cm(2),eig(2),box(2,2),work(100),detAp(maxlob)
270       double precision eig_limit /1.0D-8/
271       double precision Big /10.0D0/
272       double precision vec(3,3)
273       logical lprint,fail,lcheck
274       lcheck=.false.
275       lprint=.false.
276       fail=.false.
277       if (the.eq.0.0D0 .or. the.eq.pi) then
278 #ifdef MPI
279         write (*,'(a,i4,a,i3,a,1pe14.5)') 
280      & 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the
281 #else
282 cd        write (iout,'(a,i3,a,1pe14.5)') 
283 cd     &   'Error in GenSide: it=',it,' theta=',the
284 #endif
285         fail=.true.
286         return
287       endif
288       tant=dtan(the-pipol)
289       nlobit=nlob(it)
290       if (lprint) then
291 #ifdef MPI
292         print '(a,i4,a)','CG Processor:',me,' Enter Gen_Side.'
293         write (iout,'(a,i4,a)') 'Processor:',me,' Enter Gen_Side.'
294 #endif
295         print *,'it=',it,' nlobit=',nlobit,' the=',the,' tant=',tant
296         write (iout,*) 'it=',it,' nlobit=',nlobit,' the=',the,
297      &     ' tant=',tant
298       endif
299       do i=1,nlobit
300        zz1=tant-censc(1,i,it)
301         do k=1,3
302           do l=1,3
303             a(k,l)=gaussc(k,l,i,it)
304           enddo
305         enddo
306         detApi=a(2,2)*a(3,3)-a(2,3)**2
307         Ap_inv(2,2)=a(3,3)/detApi
308         Ap_inv(2,3)=-a(2,3)/detApi
309         Ap_inv(3,2)=Ap_inv(2,3)
310         Ap_inv(3,3)=a(2,2)/detApi
311         if (lprint) then
312           write (*,'(/a,i2/)') 'Cluster #',i
313           write (*,'(3(1pe14.5),5x,1pe14.5)') 
314      &    ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
315           write (iout,'(/a,i2/)') 'Cluster #',i
316           write (iout,'(3(1pe14.5),5x,1pe14.5)') 
317      &    ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
318         endif
319         W1i=0.0D0
320         do k=2,3
321           do l=2,3
322             W1i=W1i+a(k,1)*a(l,1)*Ap_inv(k,l)
323           enddo
324         enddo
325         W1i=a(1,1)-W1i
326         W1(i)=dexp(bsc(i,it)-0.5D0*W1i*zz1*zz1)
327 c        if (lprint) write(*,'(a,3(1pe15.5)/)')
328 c     &          'detAp, W1, anormi',detApi,W1i,anormi
329         do k=2,3
330           zk=censc(k,i,it)
331           do l=2,3
332             zk=zk+zz1*Ap_inv(k,l)*a(l,1)
333           enddo
334           z(k,i)=zk
335         enddo
336         detAp(i)=dsqrt(detApi)
337       enddo
338
339       if (lprint) then
340         print *,'W1:',(w1(i),i=1,nlobit)
341         print *,'detAp:',(detAp(i),i=1,nlobit)
342         print *,'Z'
343         do i=1,nlobit
344           print '(i2,3f10.5)',i,(rad2deg*z(j,i),j=2,3)
345         enddo
346         write (iout,*) 'W1:',(w1(i),i=1,nlobit)
347         write (iout,*) 'detAp:',(detAp(i),i=1,nlobit)
348         write (iout,*) 'Z'
349         do i=1,nlobit
350           write (iout,'(i2,3f10.5)') i,(rad2deg*z(j,i),j=2,3)
351         enddo
352       endif
353       if (lcheck) then
354 C Writing the distribution just to check the procedure
355       fac=0.0D0
356       dV=deg2rad**2*10.0D0
357       sum=0.0D0
358       sum1=0.0D0
359       do i=1,nlobit
360         fac=fac+W1(i)/detAp(i)
361       enddo 
362       fac=1.0D0/(2.0D0*fac*pi)
363 cd    print *,it,'fac=',fac
364       do ial=90,180,2
365         y(1)=deg2rad*ial
366         do iom=-180,180,5
367           y(2)=deg2rad*iom
368           wart=0.0D0
369           do i=1,nlobit
370             do j=2,3
371               do k=2,3
372                 a(j-1,k-1)=gaussc(j,k,i,it)
373               enddo
374             enddo
375             y2=y(2)
376
377             do iii=-1,1
378           
379               y(2)=y2+iii*dwapi
380
381               wykl=0.0D0
382               do j=1,2
383                 do k=1,2 
384                   wykl=wykl+a(j,k)*(y(j)-z(j+1,i))*(y(k)-z(k+1,i))
385                 enddo
386               enddo
387               wart=wart+W1(i)*dexp(-0.5D0*wykl)
388
389             enddo
390
391             y(2)=y2
392
393           enddo
394 c         print *,'y',y(1),y(2),' fac=',fac
395           wart=fac*wart
396           write (20,'(2f10.3,1pd15.5)') y(1)*rad2deg,y(2)*rad2deg,wart
397           sum=sum+wart
398           sum1=sum1+1.0D0
399         enddo
400       enddo
401 c     print *,'it=',it,' sum=',sum*dV,' sum1=',sum1*dV
402       return
403       endif
404
405 C Calculate the CM of the system
406 C
407       do i=1,nlobit
408         W1(i)=W1(i)/detAp(i)
409       enddo
410       sumW(0)=0.0D0
411       do i=1,nlobit
412         sumW(i)=sumW(i-1)+W1(i)
413       enddo
414       cm(1)=z(2,1)*W1(1)
415       cm(2)=z(3,1)*W1(1)
416       do j=2,nlobit
417         cm(1)=cm(1)+z(2,j)*W1(j) 
418         cm(2)=cm(2)+W1(j)*(z(3,1)+pinorm(z(3,j)-z(3,1)))
419       enddo
420       cm(1)=cm(1)/sumW(nlobit)
421       cm(2)=cm(2)/sumW(nlobit)
422       if (cm(1).gt.Big .or. cm(1).lt.-Big .or.
423      & cm(2).gt.Big .or. cm(2).lt.-Big) then
424 cd        write (iout,'(a)') 
425 cd     & 'Unexpected error in GenSide - CM coordinates too large.'
426 cd        write (iout,'(i5,2(1pe14.5))') it,cm(1),cm(2)
427 cd        write (*,'(a)') 
428 cd     & 'Unexpected error in GenSide - CM coordinates too large.'
429 cd        write (*,'(i5,2(1pe14.5))') it,cm(1),cm(2)
430         fail=.true. 
431         return
432       endif
433 cd    print *,'CM:',cm(1),cm(2)
434 C
435 C Find the largest search distance from CM
436 C
437       radmax=0.0D0
438       do i=1,nlobit
439         do j=2,3
440           do k=2,3
441             a(j-1,k-1)=gaussc(j,k,i,it) 
442           enddo
443         enddo
444 #ifdef NAG
445         call f02faf('N','U',2,a,3,eig,work,100,ifail)
446 #else
447         call djacob(2,3,10000,1.0d-10,a,vec,eig)
448 #endif
449 #ifdef MPI
450         if (lprint) then
451           print *,'*************** CG Processor',me
452           print *,'CM:',cm(1),cm(2)
453           write (iout,*) '*************** CG Processor',me
454           write (iout,*) 'CM:',cm(1),cm(2)
455           print '(A,8f10.5)','Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
456           write (iout,'(A,8f10.5)')
457      &        'Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
458         endif
459 #endif
460         if (eig(1).lt.eig_limit) then
461           write(iout,'(a)')
462      &     'From Mult_Norm: Eigenvalues of A are too small.'
463           write(*,'(a)')
464      &     'From Mult_Norm: Eigenvalues of A are too small.'
465           fail=.true.
466           return
467         endif
468         radius=0.0D0
469 cd      print *,'i=',i
470         do j=1,2
471           radius=radius+pinorm(z(j+1,i)-cm(j))**2
472         enddo
473         radius=dsqrt(radius)+3.0D0/dsqrt(eig(1))
474         if (radius.gt.radmax) radmax=radius
475       enddo
476       if (radmax.gt.pi) radmax=pi
477 C
478 C Determine the boundaries of the search rectangle.
479 C
480       if (lprint) then
481         print '(a,4(1pe14.4))','W1: ',(W1(i),i=1,nlob(it) )
482         print '(a,4(1pe14.4))','radmax: ',radmax
483       endif
484       box(1,1)=dmax1(cm(1)-radmax,0.0D0)
485       box(2,1)=dmin1(cm(1)+radmax,pi)
486       box(1,2)=cm(2)-radmax
487       box(2,2)=cm(2)+radmax
488       if (lprint) then
489 #ifdef MPI
490         print *,'CG Processor',me,' Array BOX:'
491 #else
492         print *,'Array BOX:'
493 #endif
494         print '(4(1pe14.4))',((box(k,j),k=1,2),j=1,2)
495         print '(a,4(1pe14.4))','sumW: ',(sumW(i),i=0,nlob(it) )
496 #ifdef MPI
497         write (iout,*)'CG Processor',me,' Array BOX:'
498 #else
499         write (iout,*)'Array BOX:'
500 #endif
501         write(iout,'(4(1pe14.4))') ((box(k,j),k=1,2),j=1,2)
502         write(iout,'(a,4(1pe14.4))')'sumW: ',(sumW(i),i=0,nlob(it) )
503       endif
504       if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
505 #ifdef MPI
506         write (iout,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
507         write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
508 #else
509 c        write (iout,'(a)') 'Bad sampling box.'
510 #endif
511         fail=.true.
512         return
513       endif
514       which_lobe=ran_number(0.0D0,sumW(nlobit))
515 c     print '(a,1pe14.4)','which_lobe=',which_lobe
516       do i=1,nlobit
517         if (sumW(i-1).le.which_lobe .and. sumW(i).ge.which_lobe) goto 1
518       enddo
519     1 ilob=i
520 c     print *,'ilob=',ilob,' nlob=',nlob(it)
521       do i=2,3
522         cm(i-1)=z(i,ilob)
523         do j=2,3
524           a(i-1,j-1)=gaussc(i,j,ilob,it)
525         enddo
526       enddo
527 cd    print '(a,i4,a)','CG Processor',me,' Calling MultNorm1.'
528       call mult_norm1(3,2,a,cm,box,y,fail)
529       if (fail) return
530       al=y(1)
531       om=pinorm(y(2))
532 cd    print *,'al=',al,' om=',om
533 cd    stop
534       return
535       end
536 c---------------------------------------------------------------------------
537       double precision function ran_number(x1,x2)
538 C Calculate a random real number from the range (x1,x2).
539       implicit real*8 (a-h,o-z)
540       include 'DIMENSIONS'
541       double precision x1,x2,fctor
542       data fctor /2147483647.0D0/
543 #ifdef MPI
544       include "mpif.h"
545       include 'COMMON.SETUP'
546       ran_number=x1+(x2-x1)*prng_next(me)
547 #else
548       call vrnd(ix,1)
549       ran_number=x1+(x2-x1)*ix/fctor
550 #endif
551       return
552       end
553 c--------------------------------------------------------------------------
554       integer function iran_num(n1,n2)
555 C Calculate a random integer number from the range (n1,n2).
556       implicit real*8 (a-h,o-z)
557       include 'DIMENSIONS'
558       integer n1,n2,ix
559       real fctor /2147483647.0/
560 #ifdef MPI
561       include "mpif.h"
562       include 'COMMON.SETUP'
563       ix=n1+(n2-n1+1)*prng_next(me)
564       if (ix.lt.n1) ix=n1
565       if (ix.gt.n2) ix=n2
566       iran_num=ix
567 #else
568       call vrnd(ix,1)
569       ix=n1+(n2-n1+1)*(ix/fctor)
570       if (ix.gt.n2) ix=n2
571       iran_num=ix
572 #endif
573       return
574       end
575 c--------------------------------------------------------------------------
576       double precision function binorm(x1,x2,sigma1,sigma2,ak)
577       implicit real*8 (a-h,o-z)
578 c     print '(a)','Enter BINORM.'
579       alowb=dmin1(x1-3.0D0*sigma1,x2-3.0D0*sigma2)
580       aupb=dmax1(x1+3.0D0*sigma1,x2+3.0D0*sigma2)
581       seg=sigma1/(sigma1+ak*sigma2)
582       alen=ran_number(0.0D0,1.0D0)
583       if (alen.lt.seg) then
584         binorm=anorm_distr(x1,sigma1,alowb,aupb)
585       else
586         binorm=anorm_distr(x2,sigma2,alowb,aupb)
587       endif
588 c     print '(a)','Exiting BINORM.'
589       return
590       end
591 c-----------------------------------------------------------------------
592 c      double precision function anorm_distr(x,sigma,alowb,aupb)
593 c      implicit real*8 (a-h,o-z)
594 c     print '(a)','Enter ANORM_DISTR.'
595 c   10 y=ran_number(alowb,aupb)
596 c      expon=dexp(-0.5D0*((y-x)/sigma)**2)
597 c      ran=ran_number(0.0D0,1.0D0)
598 c      if (expon.lt.ran) goto 10
599 c      anorm_distr=y
600 c     print '(a)','Exiting ANORM_DISTR.'
601 c      return
602 c      end
603 c-----------------------------------------------------------------------
604         double precision function anorm_distr(x,sigma,alowb,aupb)
605         implicit real*8 (a-h,o-z)
606 c  to make a normally distributed deviate with zero mean and unit variance
607 c
608         integer iset
609         real fac,gset,rsq,v1,v2,ran1
610         save iset,gset
611         data iset/0/
612         if(iset.eq.0) then
613 1               v1=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
614                 v2=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
615                 rsq=v1**2+v2**2
616                 if(rsq.ge.1.d0.or.rsq.eq.0.0d0) goto 1
617                 fac=sqrt(-2.0d0*log(rsq)/rsq)
618                 gset=v1*fac
619                 gaussdev=v2*fac
620                 iset=1
621         else
622                 gaussdev=gset
623                 iset=0
624         endif
625         anorm_distr=x+gaussdev*sigma
626       return
627       end
628 c------------------------------------------------------------------------ 
629       subroutine mult_norm(lda,n,a,x,fail)
630 C
631 C Generate the vector X whose elements obey the multiple-normal distribution
632 C from exp(-0.5*X'AX). LDA is the leading dimension of the moment matrix A,
633 C n is the dimension of the problem. FAIL is set at .TRUE., if the smallest
634 C eigenvalue of the matrix A is close to 0.
635 C
636       implicit double precision (a-h,o-z)
637       double precision a(lda,n),x(n),eig(100),vec(3,3),work(100)
638       double precision eig_limit /1.0D-8/
639       logical fail
640       fail=.false.
641 c     print '(a)','Enter MULT_NORM.'
642 C
643 C Find the smallest eigenvalue of the matrix A.
644 C
645 c     do i=1,n
646 c       print '(8f10.5)',(a(i,j),j=1,n)
647 c     enddo
648 #ifdef NAG
649       call f02faf('V','U',2,a,lda,eig,work,100,ifail)
650 #else
651       call djacob(2,lda,10000,1.0d-10,a,vec,eig)
652 #endif
653 c     print '(8f10.5)',(eig(i),i=1,n)
654 C     print '(a)'
655 c     do i=1,n
656 c       print '(8f10.5)',(a(i,j),j=1,n)
657 c     enddo
658       if (eig(1).lt.eig_limit) then
659         print *,'From Mult_Norm: Eigenvalues of A are too small.'
660         fail=.true.     
661         return
662       endif
663
664 C Generate points following the normal distributions along the principal
665 C axes of the moment matrix. Store in WORK.
666 C
667       do i=1,n
668         sigma=1.0D0/dsqrt(eig(i))
669         alim=-3.0D0*sigma
670         work(i)=anorm_distr(0.0D0,sigma,-alim,alim)
671       enddo
672 C
673 C Transform the vector of normal variables back to the original basis.
674 C
675       do i=1,n
676         xi=0.0D0
677         do j=1,n
678           xi=xi+a(i,j)*work(j)
679         enddo
680         x(i)=xi
681       enddo
682       return
683       end
684 c------------------------------------------------------------------------ 
685       subroutine mult_norm1(lda,n,a,z,box,x,fail)
686 C
687 C Generate the vector X whose elements obey the multi-gaussian multi-dimensional
688 C distribution from sum_{i=1}^m W(i)exp[-0.5*X'(i)A(i)X(i)]. LDA is the 
689 C leading dimension of the moment matrix A, n is the dimension of the 
690 C distribution, nlob is the number of lobes. FAIL is set at .TRUE., if the 
691 C smallest eigenvalue of the matrix A is close to 0.
692 C
693       implicit real*8 (a-h,o-z)
694       include 'DIMENSIONS'
695 #ifdef MPI
696       include 'mpif.h'
697 #endif
698       double precision a(lda,n),z(n),x(n),box(n,n)
699       double precision etmp
700       include 'COMMON.IOUNITS'
701 #ifdef MP
702       include 'COMMON.SETUP' 
703 #endif
704       logical fail
705
706 C Generate points following the normal distributions along the principal
707 C axes of the moment matrix. Store in WORK.
708 C
709 cd    print *,'CG Processor',me,' entered MultNorm1.'
710 cd    print '(2(1pe14.4),3x,1pe14.4)',((a(i,j),j=1,2),z(i),i=1,2)
711 cd    do i=1,n
712 cd      print *,i,box(1,i),box(2,i)
713 cd    enddo
714       istep = 0
715    10 istep = istep + 1
716       if (istep.gt.10000) then
717 c        write (iout,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
718 c     & ' in MultNorm1.'
719 c        write (*,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
720 c     & ' in MultNorm1.'
721 c        write (iout,*) 'box',box
722 c        write (iout,*) 'a',a
723 c        write (iout,*) 'z',z
724         fail=.true.
725         return
726       endif
727       do i=1,n
728        x(i)=ran_number(box(1,i),box(2,i))
729       enddo
730       ww=0.0D0
731       do i=1,n
732         xi=pinorm(x(i)-z(i))
733         ww=ww+0.5D0*a(i,i)*xi*xi
734         do j=i+1,n
735           ww=ww+a(i,j)*xi*pinorm(x(j)-z(j))
736         enddo
737       enddo
738       dec=ran_number(0.0D0,1.0D0)
739 c      print *,(x(i),i=1,n),ww,dexp(-ww),dec
740 crc   if (dec.gt.dexp(-ww)) goto 10
741       if(-ww.lt.100) then
742        etmp=dexp(-ww)
743       else
744        return  
745       endif
746       if (dec.gt.etmp) goto 10
747 cd    print *,'CG Processor',me,' exitting MultNorm1.'
748       return
749       end
750 c
751 crc--------------------------------------
752       subroutine overlap_sc(scfail)
753 c     Internal and cartesian coordinates must be consistent as input,
754 c     and will be up-to-date on return.
755 c     At the end of this procedure, scfail is true if there are
756 c     overlapping residues left, or false otherwise (success)
757       implicit real*8 (a-h,o-z)
758       include 'DIMENSIONS'
759       include 'COMMON.CHAIN'
760       include 'COMMON.INTERACT'
761       include 'COMMON.FFIELD'
762       include 'COMMON.VAR'
763       include 'COMMON.SBRIDGE'
764       include 'COMMON.IOUNITS'
765       logical had_overlaps,fail,scfail
766       integer ioverlap(maxres),ioverlap_last
767
768       had_overlaps=.false.
769       call overlap_sc_list(ioverlap,ioverlap_last)
770       if (ioverlap_last.gt.0) then
771         write (iout,*) '#OVERLAPing residues ',ioverlap_last
772         write (iout,'(20i4)') (ioverlap(k),k=1,ioverlap_last)
773         had_overlaps=.true.
774       endif
775
776       maxsi=1000
777       do k=1,1000
778         if (ioverlap_last.eq.0) exit
779
780         do ires=1,ioverlap_last 
781           i=ioverlap(ires)
782           iti=itype(i)
783           if (iti.ne.10) then
784             nsi=0
785             fail=.true.
786             do while (fail.and.nsi.le.maxsi)
787               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
788               nsi=nsi+1
789             enddo
790             if(fail) goto 999
791           endif
792         enddo
793
794         call chainbuild
795         call overlap_sc_list(ioverlap,ioverlap_last)
796 c        write (iout,*) 'Overlaping residues ',ioverlap_last,
797 c     &           (ioverlap(j),j=1,ioverlap_last)
798       enddo
799
800       if (k.le.1000.and.ioverlap_last.eq.0) then
801         scfail=.false.
802         if (had_overlaps) then
803           write (iout,*) '#OVERLAPing all corrected after ',k,
804      &         ' random generation'
805         endif
806       else
807         scfail=.true.
808         write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
809         write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last)
810       endif
811
812       return
813
814  999  continue
815       write (iout,'(a30,i5,a12,i4)') 
816      &               '#OVERLAP FAIL in gen_side after',maxsi,
817      &               'iter for RES',i
818       scfail=.true.
819       return
820       end
821
822       subroutine overlap_sc_list(ioverlap,ioverlap_last)
823       implicit real*8 (a-h,o-z)
824       include 'DIMENSIONS'
825       include 'COMMON.GEO'
826       include 'COMMON.LOCAL'
827       include 'COMMON.IOUNITS'
828       include 'COMMON.CHAIN'
829       include 'COMMON.INTERACT'
830       include 'COMMON.FFIELD'
831       include 'COMMON.VAR'
832       include 'COMMON.CALC'
833       logical fail
834       integer ioverlap(maxres),ioverlap_last
835       data redfac /0.5D0/
836
837       ioverlap_last=0
838 C Check for SC-SC overlaps and mark residues
839 c      print *,'>>overlap_sc nnt=',nnt,' nct=',nct
840       ind=0
841       do i=iatsc_s,iatsc_e
842         itypi=itype(i)
843         itypi1=itype(i+1)
844         xi=c(1,nres+i)
845         yi=c(2,nres+i)
846         zi=c(3,nres+i)
847         dxi=dc_norm(1,nres+i)
848         dyi=dc_norm(2,nres+i)
849         dzi=dc_norm(3,nres+i)
850         dsci_inv=dsc_inv(itypi)
851 c
852        do iint=1,nint_gr(i)
853          do j=istart(i,iint),iend(i,iint)
854             ind=ind+1
855             itypj=itype(j)
856             dscj_inv=dsc_inv(itypj)
857             sig0ij=sigma(itypi,itypj)
858             chi1=chi(itypi,itypj)
859             chi2=chi(itypj,itypi)
860             chi12=chi1*chi2
861             chip1=chip(itypi)
862             chip2=chip(itypj)
863             chip12=chip1*chip2
864             alf1=alp(itypi)   
865             alf2=alp(itypj)   
866             alf12=0.5D0*(alf1+alf2)
867           if (j.gt.i+1) then
868            rcomp=sigmaii(itypi,itypj)
869           else 
870            rcomp=sigma(itypi,itypj)
871           endif
872 c         print '(2(a3,2i3),a3,2f10.5)',
873 c     &        ' i=',i,iti,' j=',j,itj,' d=',dist(nres+i,nres+j)
874 c     &        ,rcomp
875             xj=c(1,nres+j)-xi
876             yj=c(2,nres+j)-yi
877             zj=c(3,nres+j)-zi
878             dxj=dc_norm(1,nres+j)
879             dyj=dc_norm(2,nres+j)
880             dzj=dc_norm(3,nres+j)
881             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
882             rij=dsqrt(rrij)
883             call sc_angular
884             sigsq=1.0D0/sigsq
885             sig=sig0ij*dsqrt(sigsq)
886             rij_shift=1.0D0/rij-sig+sig0ij
887
888 ct          if ( 1.0/rij .lt. redfac*rcomp .or. 
889 ct     &       rij_shift.le.0.0D0 ) then
890             if ( rij_shift.le.0.0D0 ) then
891 cd           write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)')
892 cd     &     'overlap SC-SC: i=',i,' j=',j,
893 cd     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
894 cd     &     rcomp,1.0/rij,rij_shift
895           ioverlap_last=ioverlap_last+1
896           ioverlap(ioverlap_last)=i         
897           do k=1,ioverlap_last-1
898            if (ioverlap(k).eq.i) ioverlap_last=ioverlap_last-1
899           enddo
900           ioverlap_last=ioverlap_last+1
901           ioverlap(ioverlap_last)=j         
902           do k=1,ioverlap_last-1
903            if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1
904           enddo 
905          endif
906         enddo
907        enddo
908       enddo
909       return
910       end