D i L aminokwasy dzialajace dla src_MD
[unres.git] / source / unres / src_MD / 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=iabs(itype(2))
19         phi(4)=gen_phi(4,iabs(itype(2)),abs(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=abs(itype(i-1))
58         it2=abs(itype(i-2))
59         it=abs(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 c         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=abs(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=abs(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=abs(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,1,1)*y(k)
242      &                +bthet(k,it,1,1)*z(k)
243       enddo
244       sig=polthet(3,it)
245       do j=2,0,-1
246         sig=sig*thet_pred_mean+polthet(j,it)
247       enddo
248       sig=0.5D0/(sig*sig+sigc0(it))
249       ak=dexp(gthet(1,it)-
250      &0.5D0*((gthet(2,it)-thet_pred_mean)/gthet(3,it))**2)
251 c     print '(i5,5(1pe14.4))',it,(gthet(j,it),j=1,3)
252 c     print '(5(1pe14.4))',thet_pred_mean,theta0(it),sig,sig0(it),ak
253       theta_temp=binorm(thet_pred_mean,theta0(it),sig,sig0(it),ak) 
254       if (theta_temp.lt.theta_min) theta_temp=theta_min
255       if (theta_temp.gt.theta_max) theta_temp=theta_max
256       gen_theta=theta_temp
257 c     print '(a)','Exiting GENTHETA.'
258       return
259       end
260 c-------------------------------------------------------------------------
261       subroutine gen_side(it,the,al,om,fail)
262       implicit real*8 (a-h,o-z)
263       include 'DIMENSIONS'
264       include 'COMMON.GEO'
265       include 'COMMON.LOCAL'
266       include 'COMMON.SETUP'
267       include 'COMMON.IOUNITS'
268       double precision MaxBoxLen /10.0D0/
269       double precision Ap_inv(3,3),a(3,3),z(3,maxlob),W1(maxlob),
270      & sumW(0:maxlob),y(2),cm(2),eig(2),box(2,2),work(100),detAp(maxlob)
271       double precision eig_limit /1.0D-8/
272       double precision Big /10.0D0/
273       double precision vec(3,3)
274       logical lprint,fail,lcheck
275       lcheck=.false.
276       lprint=.false.
277       fail=.false.
278       if (the.eq.0.0D0 .or. the.eq.pi) then
279 #ifdef MPI
280         write (*,'(a,i4,a,i3,a,1pe14.5)') 
281      & 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the
282 #else
283 cd        write (iout,'(a,i3,a,1pe14.5)') 
284 cd     &   'Error in GenSide: it=',it,' theta=',the
285 #endif
286         fail=.true.
287         return
288       endif
289       tant=dtan(the-pipol)
290       nlobit=nlob(it)
291       if (lprint) then
292 #ifdef MPI
293         print '(a,i4,a)','CG Processor:',me,' Enter Gen_Side.'
294         write (iout,'(a,i4,a)') 'Processor:',me,' Enter Gen_Side.'
295 #endif
296         print *,'it=',it,' nlobit=',nlobit,' the=',the,' tant=',tant
297         write (iout,*) 'it=',it,' nlobit=',nlobit,' the=',the,
298      &     ' tant=',tant
299       endif
300       do i=1,nlobit
301        zz1=tant-censc(1,i,it)
302         do k=1,3
303           do l=1,3
304             a(k,l)=gaussc(k,l,i,it)
305           enddo
306         enddo
307         detApi=a(2,2)*a(3,3)-a(2,3)**2
308         Ap_inv(2,2)=a(3,3)/detApi
309         Ap_inv(2,3)=-a(2,3)/detApi
310         Ap_inv(3,2)=Ap_inv(2,3)
311         Ap_inv(3,3)=a(2,2)/detApi
312         if (lprint) then
313           write (*,'(/a,i2/)') 'Cluster #',i
314           write (*,'(3(1pe14.5),5x,1pe14.5)') 
315      &    ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
316           write (iout,'(/a,i2/)') 'Cluster #',i
317           write (iout,'(3(1pe14.5),5x,1pe14.5)') 
318      &    ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
319         endif
320         W1i=0.0D0
321         do k=2,3
322           do l=2,3
323             W1i=W1i+a(k,1)*a(l,1)*Ap_inv(k,l)
324           enddo
325         enddo
326         W1i=a(1,1)-W1i
327         W1(i)=dexp(bsc(i,it)-0.5D0*W1i*zz1*zz1)
328 c        if (lprint) write(*,'(a,3(1pe15.5)/)')
329 c     &          'detAp, W1, anormi',detApi,W1i,anormi
330         do k=2,3
331           zk=censc(k,i,it)
332           do l=2,3
333             zk=zk+zz1*Ap_inv(k,l)*a(l,1)
334           enddo
335           z(k,i)=zk
336         enddo
337         detAp(i)=dsqrt(detApi)
338       enddo
339
340       if (lprint) then
341         print *,'W1:',(w1(i),i=1,nlobit)
342         print *,'detAp:',(detAp(i),i=1,nlobit)
343         print *,'Z'
344         do i=1,nlobit
345           print '(i2,3f10.5)',i,(rad2deg*z(j,i),j=2,3)
346         enddo
347         write (iout,*) 'W1:',(w1(i),i=1,nlobit)
348         write (iout,*) 'detAp:',(detAp(i),i=1,nlobit)
349         write (iout,*) 'Z'
350         do i=1,nlobit
351           write (iout,'(i2,3f10.5)') i,(rad2deg*z(j,i),j=2,3)
352         enddo
353       endif
354       if (lcheck) then
355 C Writing the distribution just to check the procedure
356       fac=0.0D0
357       dV=deg2rad**2*10.0D0
358       sum=0.0D0
359       sum1=0.0D0
360       do i=1,nlobit
361         fac=fac+W1(i)/detAp(i)
362       enddo 
363       fac=1.0D0/(2.0D0*fac*pi)
364 cd    print *,it,'fac=',fac
365       do ial=90,180,2
366         y(1)=deg2rad*ial
367         do iom=-180,180,5
368           y(2)=deg2rad*iom
369           wart=0.0D0
370           do i=1,nlobit
371             do j=2,3
372               do k=2,3
373                 a(j-1,k-1)=gaussc(j,k,i,it)
374               enddo
375             enddo
376             y2=y(2)
377
378             do iii=-1,1
379           
380               y(2)=y2+iii*dwapi
381
382               wykl=0.0D0
383               do j=1,2
384                 do k=1,2 
385                   wykl=wykl+a(j,k)*(y(j)-z(j+1,i))*(y(k)-z(k+1,i))
386                 enddo
387               enddo
388               wart=wart+W1(i)*dexp(-0.5D0*wykl)
389
390             enddo
391
392             y(2)=y2
393
394           enddo
395 c         print *,'y',y(1),y(2),' fac=',fac
396           wart=fac*wart
397           write (20,'(2f10.3,1pd15.5)') y(1)*rad2deg,y(2)*rad2deg,wart
398           sum=sum+wart
399           sum1=sum1+1.0D0
400         enddo
401       enddo
402 c     print *,'it=',it,' sum=',sum*dV,' sum1=',sum1*dV
403       return
404       endif
405
406 C Calculate the CM of the system
407 C
408       do i=1,nlobit
409         W1(i)=W1(i)/detAp(i)
410       enddo
411       sumW(0)=0.0D0
412       do i=1,nlobit
413         sumW(i)=sumW(i-1)+W1(i)
414       enddo
415       cm(1)=z(2,1)*W1(1)
416       cm(2)=z(3,1)*W1(1)
417       do j=2,nlobit
418         cm(1)=cm(1)+z(2,j)*W1(j) 
419         cm(2)=cm(2)+W1(j)*(z(3,1)+pinorm(z(3,j)-z(3,1)))
420       enddo
421       cm(1)=cm(1)/sumW(nlobit)
422       cm(2)=cm(2)/sumW(nlobit)
423       if (cm(1).gt.Big .or. cm(1).lt.-Big .or.
424      & cm(2).gt.Big .or. cm(2).lt.-Big) then
425 cd        write (iout,'(a)') 
426 cd     & 'Unexpected error in GenSide - CM coordinates too large.'
427 cd        write (iout,'(i5,2(1pe14.5))') it,cm(1),cm(2)
428 cd        write (*,'(a)') 
429 cd     & 'Unexpected error in GenSide - CM coordinates too large.'
430 cd        write (*,'(i5,2(1pe14.5))') it,cm(1),cm(2)
431         fail=.true. 
432         return
433       endif
434 cd    print *,'CM:',cm(1),cm(2)
435 C
436 C Find the largest search distance from CM
437 C
438       radmax=0.0D0
439       do i=1,nlobit
440         do j=2,3
441           do k=2,3
442             a(j-1,k-1)=gaussc(j,k,i,it) 
443           enddo
444         enddo
445 #ifdef NAG
446         call f02faf('N','U',2,a,3,eig,work,100,ifail)
447 #else
448         call djacob(2,3,10000,1.0d-10,a,vec,eig)
449 #endif
450 #ifdef MPI
451         if (lprint) then
452           print *,'*************** CG Processor',me
453           print *,'CM:',cm(1),cm(2)
454           write (iout,*) '*************** CG Processor',me
455           write (iout,*) 'CM:',cm(1),cm(2)
456           print '(A,8f10.5)','Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
457           write (iout,'(A,8f10.5)')
458      &        'Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
459         endif
460 #endif
461         if (eig(1).lt.eig_limit) then
462           write(iout,'(a)')
463      &     'From Mult_Norm: Eigenvalues of A are too small.'
464           write(*,'(a)')
465      &     'From Mult_Norm: Eigenvalues of A are too small.'
466           fail=.true.
467           return
468         endif
469         radius=0.0D0
470 cd      print *,'i=',i
471         do j=1,2
472           radius=radius+pinorm(z(j+1,i)-cm(j))**2
473         enddo
474         radius=dsqrt(radius)+3.0D0/dsqrt(eig(1))
475         if (radius.gt.radmax) radmax=radius
476       enddo
477       if (radmax.gt.pi) radmax=pi
478 C
479 C Determine the boundaries of the search rectangle.
480 C
481       if (lprint) then
482         print '(a,4(1pe14.4))','W1: ',(W1(i),i=1,nlob(it) )
483         print '(a,4(1pe14.4))','radmax: ',radmax
484       endif
485       box(1,1)=dmax1(cm(1)-radmax,0.0D0)
486       box(2,1)=dmin1(cm(1)+radmax,pi)
487       box(1,2)=cm(2)-radmax
488       box(2,2)=cm(2)+radmax
489       if (lprint) then
490 #ifdef MPI
491         print *,'CG Processor',me,' Array BOX:'
492 #else
493         print *,'Array BOX:'
494 #endif
495         print '(4(1pe14.4))',((box(k,j),k=1,2),j=1,2)
496         print '(a,4(1pe14.4))','sumW: ',(sumW(i),i=0,nlob(it) )
497 #ifdef MPI
498         write (iout,*)'CG Processor',me,' Array BOX:'
499 #else
500         write (iout,*)'Array BOX:'
501 #endif
502         write(iout,'(4(1pe14.4))') ((box(k,j),k=1,2),j=1,2)
503         write(iout,'(a,4(1pe14.4))')'sumW: ',(sumW(i),i=0,nlob(it) )
504       endif
505       if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
506 #ifdef MPI
507         write (iout,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
508         write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
509 #else
510 c        write (iout,'(a)') 'Bad sampling box.'
511 #endif
512         fail=.true.
513         return
514       endif
515       which_lobe=ran_number(0.0D0,sumW(nlobit))
516 c     print '(a,1pe14.4)','which_lobe=',which_lobe
517       do i=1,nlobit
518         if (sumW(i-1).le.which_lobe .and. sumW(i).ge.which_lobe) goto 1
519       enddo
520     1 ilob=i
521 c     print *,'ilob=',ilob,' nlob=',nlob(it)
522       do i=2,3
523         cm(i-1)=z(i,ilob)
524         do j=2,3
525           a(i-1,j-1)=gaussc(i,j,ilob,it)
526         enddo
527       enddo
528 cd    print '(a,i4,a)','CG Processor',me,' Calling MultNorm1.'
529       call mult_norm1(3,2,a,cm,box,y,fail)
530       if (fail) return
531       al=y(1)
532       om=pinorm(y(2))
533 cd    print *,'al=',al,' om=',om
534 cd    stop
535       return
536       end
537 c---------------------------------------------------------------------------
538       double precision function ran_number(x1,x2)
539 C Calculate a random real number from the range (x1,x2).
540       implicit real*8 (a-h,o-z)
541       include 'DIMENSIONS'
542       double precision x1,x2,fctor
543       data fctor /2147483647.0D0/
544 #ifdef MPI
545       include "mpif.h"
546       include 'COMMON.SETUP'
547       ran_number=x1+(x2-x1)*prng_next(me)
548 #else
549       call vrnd(ix,1)
550       ran_number=x1+(x2-x1)*ix/fctor
551 #endif
552       return
553       end
554 c--------------------------------------------------------------------------
555       integer function iran_num(n1,n2)
556 C Calculate a random integer number from the range (n1,n2).
557       implicit real*8 (a-h,o-z)
558       include 'DIMENSIONS'
559       integer n1,n2,ix
560       real fctor /2147483647.0/
561 #ifdef MPI
562       include "mpif.h"
563       include 'COMMON.SETUP'
564       ix=n1+(n2-n1+1)*prng_next(me)
565       if (ix.lt.n1) ix=n1
566       if (ix.gt.n2) ix=n2
567       iran_num=ix
568 #else
569       call vrnd(ix,1)
570       ix=n1+(n2-n1+1)*(ix/fctor)
571       if (ix.gt.n2) ix=n2
572       iran_num=ix
573 #endif
574       return
575       end
576 c--------------------------------------------------------------------------
577       double precision function binorm(x1,x2,sigma1,sigma2,ak)
578       implicit real*8 (a-h,o-z)
579 c     print '(a)','Enter BINORM.'
580       alowb=dmin1(x1-3.0D0*sigma1,x2-3.0D0*sigma2)
581       aupb=dmax1(x1+3.0D0*sigma1,x2+3.0D0*sigma2)
582       seg=sigma1/(sigma1+ak*sigma2)
583       alen=ran_number(0.0D0,1.0D0)
584       if (alen.lt.seg) then
585         binorm=anorm_distr(x1,sigma1,alowb,aupb)
586       else
587         binorm=anorm_distr(x2,sigma2,alowb,aupb)
588       endif
589 c     print '(a)','Exiting BINORM.'
590       return
591       end
592 c-----------------------------------------------------------------------
593 c      double precision function anorm_distr(x,sigma,alowb,aupb)
594 c      implicit real*8 (a-h,o-z)
595 c     print '(a)','Enter ANORM_DISTR.'
596 c   10 y=ran_number(alowb,aupb)
597 c      expon=dexp(-0.5D0*((y-x)/sigma)**2)
598 c      ran=ran_number(0.0D0,1.0D0)
599 c      if (expon.lt.ran) goto 10
600 c      anorm_distr=y
601 c     print '(a)','Exiting ANORM_DISTR.'
602 c      return
603 c      end
604 c-----------------------------------------------------------------------
605         double precision function anorm_distr(x,sigma,alowb,aupb)
606         implicit real*8 (a-h,o-z)
607 c  to make a normally distributed deviate with zero mean and unit variance
608 c
609         integer iset
610         real fac,gset,rsq,v1,v2,ran1
611         save iset,gset
612         data iset/0/
613         if(iset.eq.0) then
614 1               v1=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
615                 v2=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
616                 rsq=v1**2+v2**2
617                 if(rsq.ge.1.d0.or.rsq.eq.0.0d0) goto 1
618                 fac=sqrt(-2.0d0*log(rsq)/rsq)
619                 gset=v1*fac
620                 gaussdev=v2*fac
621                 iset=1
622         else
623                 gaussdev=gset
624                 iset=0
625         endif
626         anorm_distr=x+gaussdev*sigma
627       return
628       end
629 c------------------------------------------------------------------------ 
630       subroutine mult_norm(lda,n,a,x,fail)
631 C
632 C Generate the vector X whose elements obey the multiple-normal distribution
633 C from exp(-0.5*X'AX). LDA is the leading dimension of the moment matrix A,
634 C n is the dimension of the problem. FAIL is set at .TRUE., if the smallest
635 C eigenvalue of the matrix A is close to 0.
636 C
637       implicit double precision (a-h,o-z)
638       double precision a(lda,n),x(n),eig(100),vec(3,3),work(100)
639       double precision eig_limit /1.0D-8/
640       logical fail
641       fail=.false.
642 c     print '(a)','Enter MULT_NORM.'
643 C
644 C Find the smallest eigenvalue of the matrix A.
645 C
646 c     do i=1,n
647 c       print '(8f10.5)',(a(i,j),j=1,n)
648 c     enddo
649 #ifdef NAG
650       call f02faf('V','U',2,a,lda,eig,work,100,ifail)
651 #else
652       call djacob(2,lda,10000,1.0d-10,a,vec,eig)
653 #endif
654 c     print '(8f10.5)',(eig(i),i=1,n)
655 C     print '(a)'
656 c     do i=1,n
657 c       print '(8f10.5)',(a(i,j),j=1,n)
658 c     enddo
659       if (eig(1).lt.eig_limit) then
660         print *,'From Mult_Norm: Eigenvalues of A are too small.'
661         fail=.true.     
662         return
663       endif
664
665 C Generate points following the normal distributions along the principal
666 C axes of the moment matrix. Store in WORK.
667 C
668       do i=1,n
669         sigma=1.0D0/dsqrt(eig(i))
670         alim=-3.0D0*sigma
671         work(i)=anorm_distr(0.0D0,sigma,-alim,alim)
672       enddo
673 C
674 C Transform the vector of normal variables back to the original basis.
675 C
676       do i=1,n
677         xi=0.0D0
678         do j=1,n
679           xi=xi+a(i,j)*work(j)
680         enddo
681         x(i)=xi
682       enddo
683       return
684       end
685 c------------------------------------------------------------------------ 
686       subroutine mult_norm1(lda,n,a,z,box,x,fail)
687 C
688 C Generate the vector X whose elements obey the multi-gaussian multi-dimensional
689 C distribution from sum_{i=1}^m W(i)exp[-0.5*X'(i)A(i)X(i)]. LDA is the 
690 C leading dimension of the moment matrix A, n is the dimension of the 
691 C distribution, nlob is the number of lobes. FAIL is set at .TRUE., if the 
692 C smallest eigenvalue of the matrix A is close to 0.
693 C
694       implicit real*8 (a-h,o-z)
695       include 'DIMENSIONS'
696 #ifdef MPI
697       include 'mpif.h'
698 #endif
699       double precision a(lda,n),z(n),x(n),box(n,n)
700       double precision etmp
701       include 'COMMON.IOUNITS'
702 #ifdef MP
703       include 'COMMON.SETUP' 
704 #endif
705       logical fail
706
707 C Generate points following the normal distributions along the principal
708 C axes of the moment matrix. Store in WORK.
709 C
710 cd    print *,'CG Processor',me,' entered MultNorm1.'
711 cd    print '(2(1pe14.4),3x,1pe14.4)',((a(i,j),j=1,2),z(i),i=1,2)
712 cd    do i=1,n
713 cd      print *,i,box(1,i),box(2,i)
714 cd    enddo
715       istep = 0
716    10 istep = istep + 1
717       if (istep.gt.10000) then
718 c        write (iout,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
719 c     & ' in MultNorm1.'
720 c        write (*,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
721 c     & ' in MultNorm1.'
722 c        write (iout,*) 'box',box
723 c        write (iout,*) 'a',a
724 c        write (iout,*) 'z',z
725         fail=.true.
726         return
727       endif
728       do i=1,n
729        x(i)=ran_number(box(1,i),box(2,i))
730       enddo
731       ww=0.0D0
732       do i=1,n
733         xi=pinorm(x(i)-z(i))
734         ww=ww+0.5D0*a(i,i)*xi*xi
735         do j=i+1,n
736           ww=ww+a(i,j)*xi*pinorm(x(j)-z(j))
737         enddo
738       enddo
739       dec=ran_number(0.0D0,1.0D0)
740 c      print *,(x(i),i=1,n),ww,dexp(-ww),dec
741 crc   if (dec.gt.dexp(-ww)) goto 10
742       if(-ww.lt.100) then
743        etmp=dexp(-ww)
744       else
745        return  
746       endif
747       if (dec.gt.etmp) goto 10
748 cd    print *,'CG Processor',me,' exitting MultNorm1.'
749       return
750       end
751 c
752 crc--------------------------------------
753       subroutine overlap_sc(scfail)
754 c     Internal and cartesian coordinates must be consistent as input,
755 c     and will be up-to-date on return.
756 c     At the end of this procedure, scfail is true if there are
757 c     overlapping residues left, or false otherwise (success)
758       implicit real*8 (a-h,o-z)
759       include 'DIMENSIONS'
760       include 'COMMON.CHAIN'
761       include 'COMMON.INTERACT'
762       include 'COMMON.FFIELD'
763       include 'COMMON.VAR'
764       include 'COMMON.SBRIDGE'
765       include 'COMMON.IOUNITS'
766       logical had_overlaps,fail,scfail
767       integer ioverlap(maxres),ioverlap_last
768
769       had_overlaps=.false.
770       call overlap_sc_list(ioverlap,ioverlap_last)
771       if (ioverlap_last.gt.0) then
772         write (iout,*) '#OVERLAPing residues ',ioverlap_last
773         write (iout,'(20i4)') (ioverlap(k),k=1,ioverlap_last)
774         had_overlaps=.true.
775       endif
776
777       maxsi=1000
778       do k=1,1000
779         if (ioverlap_last.eq.0) exit
780
781         do ires=1,ioverlap_last 
782           i=ioverlap(ires)
783           iti=abs(itype(i))
784           if (iti.ne.10) then
785             nsi=0
786             fail=.true.
787             do while (fail.and.nsi.le.maxsi)
788               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
789               nsi=nsi+1
790             enddo
791             if(fail) goto 999
792           endif
793         enddo
794
795         call chainbuild
796         call overlap_sc_list(ioverlap,ioverlap_last)
797 c        write (iout,*) 'Overlaping residues ',ioverlap_last,
798 c     &           (ioverlap(j),j=1,ioverlap_last)
799       enddo
800
801       if (k.le.1000.and.ioverlap_last.eq.0) then
802         scfail=.false.
803         if (had_overlaps) then
804           write (iout,*) '#OVERLAPing all corrected after ',k,
805      &         ' random generation'
806         endif
807       else
808         scfail=.true.
809         write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
810         write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last)
811       endif
812
813       return
814
815  999  continue
816       write (iout,'(a30,i5,a12,i4)') 
817      &               '#OVERLAP FAIL in gen_side after',maxsi,
818      &               'iter for RES',i
819       scfail=.true.
820       return
821       end
822
823       subroutine overlap_sc_list(ioverlap,ioverlap_last)
824       implicit real*8 (a-h,o-z)
825       include 'DIMENSIONS'
826       include 'COMMON.GEO'
827       include 'COMMON.LOCAL'
828       include 'COMMON.IOUNITS'
829       include 'COMMON.CHAIN'
830       include 'COMMON.INTERACT'
831       include 'COMMON.FFIELD'
832       include 'COMMON.VAR'
833       include 'COMMON.CALC'
834       logical fail
835       integer ioverlap(maxres),ioverlap_last
836       data redfac /0.5D0/
837
838       ioverlap_last=0
839 C Check for SC-SC overlaps and mark residues
840 c      print *,'>>overlap_sc nnt=',nnt,' nct=',nct
841       ind=0
842       do i=iatsc_s,iatsc_e
843         itypi=abs(itype(i))
844         itypi1=abs(itype(i+1))
845         xi=c(1,nres+i)
846         yi=c(2,nres+i)
847         zi=c(3,nres+i)
848         dxi=dc_norm(1,nres+i)
849         dyi=dc_norm(2,nres+i)
850         dzi=dc_norm(3,nres+i)
851         dsci_inv=dsc_inv(itypi)
852 c
853        do iint=1,nint_gr(i)
854          do j=istart(i,iint),iend(i,iint)
855             ind=ind+1
856             itypj=abs(itype(j))
857             dscj_inv=dsc_inv(itypj)
858             sig0ij=sigma(itypi,itypj)
859             chi1=chi(itypi,itypj)
860             chi2=chi(itypj,itypi)
861             chi12=chi1*chi2
862             chip1=chip(itypi)
863             chip2=chip(itypj)
864             chip12=chip1*chip2
865             alf1=alp(itypi)   
866             alf2=alp(itypj)   
867             alf12=0.5D0*(alf1+alf2)
868           if (j.gt.i+1) then
869            rcomp=sigmaii(itypi,itypj)
870           else 
871            rcomp=sigma(itypi,itypj)
872           endif
873 c         print '(2(a3,2i3),a3,2f10.5)',
874 c     &        ' i=',i,iti,' j=',j,itj,' d=',dist(nres+i,nres+j)
875 c     &        ,rcomp
876             xj=c(1,nres+j)-xi
877             yj=c(2,nres+j)-yi
878             zj=c(3,nres+j)-zi
879             dxj=dc_norm(1,nres+j)
880             dyj=dc_norm(2,nres+j)
881             dzj=dc_norm(3,nres+j)
882             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
883             rij=dsqrt(rrij)
884             call sc_angular
885             sigsq=1.0D0/sigsq
886             sig=sig0ij*dsqrt(sigsq)
887             rij_shift=1.0D0/rij-sig+sig0ij
888
889 ct          if ( 1.0/rij .lt. redfac*rcomp .or. 
890 ct     &       rij_shift.le.0.0D0 ) then
891             if ( rij_shift.le.0.0D0 ) then
892 cd           write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)')
893 cd     &     'overlap SC-SC: i=',i,' j=',j,
894 cd     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
895 cd     &     rcomp,1.0/rij,rij_shift
896           ioverlap_last=ioverlap_last+1
897           ioverlap(ioverlap_last)=i         
898           do k=1,ioverlap_last-1
899            if (ioverlap(k).eq.i) ioverlap_last=ioverlap_last-1
900           enddo
901           ioverlap_last=ioverlap_last+1
902           ioverlap(ioverlap_last)=j         
903           do k=1,ioverlap_last-1
904            if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1
905           enddo 
906          endif
907         enddo
908        enddo
909       enddo
910       return
911       end