Sept 20/2012 by Adam: commit changes mostly in the MINIM version

diff --git a/bin/unres/MD/unres_ifort_MPICH_E0LL2Y-test.exe b/bin/unres/MD/unres_ifort_MPICH_E0LL2Y-test.exe
new file mode 100755 (executable)
index 0000000..e3428c6
Binary files /dev/null and b/bin/unres/MD/unres_ifort_MPICH_E0LL2Y-test.exe differ

diff --git a/bin/unres/MD/unres_ifort_MPICH_E0LL2Y.exe b/bin/unres/MD/unres_ifort_MPICH_E0LL2Y.exe
new file mode 100755 (executable)
index 0000000..2a360a6
Binary files /dev/null and b/bin/unres/MD/unres_ifort_MPICH_E0LL2Y.exe differ

diff --git a/bin/unres/MINIM/unres_ifort_MIN_single_E0LL2Y.exe b/bin/unres/MINIM/unres_ifort_MIN_single_E0LL2Y.exe
new file mode 100755 (executable)
index 0000000..cb1a624
Binary files /dev/null and b/bin/unres/MINIM/unres_ifort_MIN_single_E0LL2Y.exe differ

diff --git a/bin/unres/MINIM/unres_ifort_MIN_single_GAB.exe b/bin/unres/MINIM/unres_ifort_MIN_single_GAB.exe
new file mode 100755 (executable)
index 0000000..394eaac
Binary files /dev/null and b/bin/unres/MINIM/unres_ifort_MIN_single_GAB.exe differ

diff --git a/bin/wham/wham_ifort_MPICH_E0LL2Y-DEBUG.exe b/bin/wham/wham_ifort_MPICH_E0LL2Y-DEBUG.exe
new file mode 100755 (executable)
index 0000000..a9af08d
Binary files /dev/null and b/bin/wham/wham_ifort_MPICH_E0LL2Y-DEBUG.exe differ

diff --git a/bin/wham/wham_ifort_MPICH_GAB-DEBUG.exe b/bin/wham/wham_ifort_MPICH_GAB-DEBUG.exe
new file mode 100755 (executable)
index 0000000..c85f077
Binary files /dev/null and b/bin/wham/wham_ifort_MPICH_GAB-DEBUG.exe differ

diff --git a/bin/wham/wham_ifort_MPICH_GAB.exe b/bin/wham/wham_ifort_MPICH_GAB.exe
new file mode 100755 (executable)
index 0000000..e1c4ae5
Binary files /dev/null and b/bin/wham/wham_ifort_MPICH_GAB.exe differ

diff --git a/source/unres/src_MIN/Makefile_gfortran_single b/source/unres/src_MIN/Makefile_gfortran_single
new file mode 100644 (file)
index 0000000..5701c47
--- /dev/null
+++ b/source/unres/src_MIN/Makefile_gfortran_single
@@ -0,0 +1,88 @@
+FC= gfortran
+
+OPT =  -O
+
+FFLAGS = -c ${OPT}
+FFLAGS1 = -c -g -C
+FFLAGS2 = -c -g -O0
+FFLAGSE = -c -O3
+
+LIBS =
+
+ARCH = LINUX
+PP = /lib/cpp -P
+
+
+all: unres
+
+.SUFFIXES: .F.f
+.F.o:
+       ${FC} ${FFLAGS} ${CPPFLAGS} $*.F
+.f.o:
+       ${FC} ${FFLAGS} ${CPPFLAGS} $*.f
+
+
+object = unres_min.o arcos.o cartprint.o chainbuild.o initialize_p.o \
+        matmult.o readrtns_min.o parmread.o \
+        pinorm.o randgens.o rescode.o intcor.o timing.o misc.o intlocal.o \
+        cartder.o checkder_p.o econstr_local.o energy_p_new_barrier.o \
+       gradient_p.o minimize_p.o sumsld.o \
+        cored.o rmdd.o geomout_min.o readpdb.o \
+        intcartderiv.o \
+        MP.o printmat.o convert.o int_to_cart.o \
+       djacob.o gen_rand_conf.o sc_move.o
+
+GAB: CPPFLAGS = -DPROCOR -DLINUX -DUNRES -DISNAN \
+       -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC
+GAB: BIN = ../../../bin/unres/MINIM/unres_gfortran_MIN_single_GAB.exe
+GAB: ${object} 
+       cc -o compinfo compinfo.c
+       ./compinfo | true
+       ${FC} ${FFLAGS} cinfo.f
+       ${FC} ${OPT} ${object} cinfo.o ${LIBS}  -o ${BIN}
+
+E0LL2Y: CPPFLAGS = -DPROCOR -DLINUX -DUNRES -DISNAN \
+       -DSPLITELE -DLANG0
+E0LL2Y: BIN = ../../../bin/unres/MINIM/unres_gfortran_MIN_single_E0LL2Y.exe
+E0LL2Y: ${object} 
+       cc -o compinfo compinfo.c
+       ./compinfo | true
+       ${FC} ${FFLAGS} cinfo.f
+       ${FC} ${OPT} ${object} cinfo.o ${LIBS}  -o ${BIN}
+
+
+clean:
+       /bin/rm *.o
+
+chainbuild.o: chainbuild.F
+       ${FC} ${FFLAGS} ${CPPFLAGS} chainbuild.F
+
+matmult.o: matmult.f
+       ${FC} ${FFLAGS} ${CPPFLAGS} matmult.f
+
+parmread.o : parmread.F
+       ${FC} ${FFLAGS} ${CPPFLAGS} parmread.F
+
+intcor.o : intcor.f
+       ${FC} ${FFLAGS} ${CPPFLAGS} intcor.f
+
+cartder.o : cartder.F
+       ${FC} ${FFLAGS} ${CPPFLAGS} cartder.F
+
+readpdb.o : readpdb.F
+       ${FC} ${FFLAGS2} ${CPPFLAGS} readpdb.F
+
+sumsld.o : sumsld.f
+       ${FC} ${FFLAGS} ${CPPFLAGS} sumsld.f
+        
+cored.o : cored.f
+       ${FC} ${FFLAGS1} ${CPPFLAGS} cored.f
+ 
+rmdd.o : rmdd.f
+       ${FC} ${FFLAGS} ${CPPFLAGS} rmdd.f
+
+energy_p_new_barrier.o : energy_p_new_barrier.F
+       ${FC} ${FFLAGSE} ${CPPFLAGS} energy_p_new_barrier.F
+
+gradient_p.o : gradient_p.F
+       ${FC} ${FFLAGSE} ${CPPFLAGS} gradient_p.F

diff --git a/source/unres/src_MIN/Makefile_ifort_single b/source/unres/src_MIN/Makefile_ifort_single
new file mode 100644 (file)
index 0000000..d21cc98
--- /dev/null
+++ b/source/unres/src_MIN/Makefile_ifort_single
@@ -0,0 +1,88 @@
+FC= ifort -g
+
+OPT =  -O3 -ip -w 
+
+FFLAGS = -c ${OPT} 
+FFLAGS1 = -c -w -g -d2 -CA -CB 
+FFLAGS2 = -c -w -g -O0 
+FFLAGSE = -c -w -O3 -ipo -ipo_obj  -opt_report 
+
+LIBS =  -lpthread 
+
+ARCH = LINUX
+PP = /lib/cpp -P
+
+
+all: unres
+
+.SUFFIXES: .F.f
+.F.o:
+       ${FC} ${FFLAGS} ${CPPFLAGS} $*.F
+.f.o:
+       ${FC} ${FFLAGS} ${CPPFLAGS} $*.f
+
+
+object = unres_min.o arcos.o cartprint.o chainbuild.o initialize_p.o \
+        matmult.o readrtns_min.o parmread.o \
+        pinorm.o randgens.o rescode.o intcor.o timing.o misc.o intlocal.o \
+        cartder.o checkder_p.o econstr_local.o energy_p_new_barrier.o \
+       gradient_p.o minimize_p.o sumsld.o \
+        cored.o rmdd.o geomout_min.o readpdb.o \
+        intcartderiv.o \
+        MP.o printmat.o convert.o int_to_cart.o \
+       djacob.o gen_rand_conf.o sc_move.o
+
+GAB: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNRES -DISNAN \
+       -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC
+GAB: BIN = ../../../bin/unres/MINIM/unres_ifort_MIN_single_GAB.exe
+GAB: ${object} 
+       cc -o compinfo compinfo.c
+       ./compinfo | true
+       ${FC} ${FFLAGS} cinfo.f
+       ${FC} ${OPT} ${object} cinfo.o ${LIBS}  -o ${BIN}
+
+E0LL2Y: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNRES -DISNAN \
+       -DSPLITELE -DLANG0
+E0LL2Y: BIN = ../../../bin/unres/MINIM/unres_ifort_MIN_single_E0LL2Y.exe
+E0LL2Y: ${object} 
+       cc -o compinfo compinfo.c
+       ./compinfo | true
+       ${FC} ${FFLAGS} cinfo.f
+       ${FC} ${OPT} ${object} cinfo.o ${LIBS}  -o ${BIN}
+
+
+clean:
+       /bin/rm *.o
+
+chainbuild.o: chainbuild.F
+       ${FC} ${FFLAGS} ${CPPFLAGS} chainbuild.F
+
+matmult.o: matmult.f
+       ${FC} ${FFLAGS} ${CPPFLAGS} matmult.f
+
+parmread.o : parmread.F
+       ${FC} ${FFLAGS} ${CPPFLAGS} parmread.F
+
+intcor.o : intcor.f
+       ${FC} ${FFLAGS} ${CPPFLAGS} intcor.f
+
+cartder.o : cartder.F
+       ${FC} ${FFLAGS} ${CPPFLAGS} cartder.F
+
+readpdb.o : readpdb.F
+       ${FC} ${FFLAGS2} ${CPPFLAGS} readpdb.F
+
+sumsld.o : sumsld.f
+       ${FC} ${FFLAGS} ${CPPFLAGS} sumsld.f
+        
+cored.o : cored.f
+       ${FC} ${FFLAGS1} ${CPPFLAGS} cored.f
+ 
+rmdd.o : rmdd.f
+       ${FC} ${FFLAGS} ${CPPFLAGS} rmdd.f
+
+energy_p_new_barrier.o : energy_p_new_barrier.F
+       ${FC} ${FFLAGSE} ${CPPFLAGS} energy_p_new_barrier.F
+
+gradient_p.o : gradient_p.F
+       ${FC} ${FFLAGSE} ${CPPFLAGS} gradient_p.F

diff --git a/source/unres/src_MIN/compinfo.c b/source/unres/src_MIN/compinfo.c
new file mode 100644 (file)
index 0000000..e28f686
--- /dev/null
+++ b/source/unres/src_MIN/compinfo.c
@@ -0,0 +1,82 @@
+#include <stdio.h>
+#include <sys/utsname.h>
+#include <sys/types.h>
+#include <time.h>
+#include <string.h>
+
+main()
+{
+FILE *in, *in1, *out;
+int i,j,k,iv1,iv2,iv3;
+char *p1,buf[500],buf1[500],buf2[100],buf3[100];
+struct utsname Name;
+time_t Tp;
+
+in=fopen("cinfo.f","r");
+out=fopen("cinfo.f.new","w");
+if (fgets(buf,498,in) != NULL)
+       fprintf(out,"C DO NOT EDIT THIS FILE - IT HAS BEEN GENERATED BY COMPINFO.C\n");
+if (fgets(buf,498,in) != NULL)
+       sscanf(&buf[1],"%d %d %d",&iv1,&iv2,&iv3);
+iv3++;
+fprintf(out,"C %d %d %d\n",iv1,iv2,iv3);
+fprintf(out,"      subroutine cinfo\n");
+fprintf(out,"      include 'COMMON.IOUNITS'\n");
+fprintf(out,"      write(iout,*)'++++ Compile info ++++'\n");
+fprintf(out,"      write(iout,*)'Version %d.%-d build %d'\n",iv1,iv2,iv3);
+uname(&Name);
+time(&Tp);
+system("whoami > tmptmp");
+in1=fopen("tmptmp","r");
+if (fscanf(in1,"%s",buf1) != EOF)
+{
+p1=ctime(&Tp);
+p1[strlen(p1)-1]='\0';
+fprintf(out,"      write(iout,*)'compiled %s'\n",p1);
+fprintf(out,"      write(iout,*)'compiled by %s@%s'\n",buf1,Name.nodename);
+fprintf(out,"      write(iout,*)'OS name:    %s '\n",Name.sysname);
+fprintf(out,"      write(iout,*)'OS release: %s '\n",Name.release);
+fprintf(out,"      write(iout,*)'OS version:',\n");
+fprintf(out,"     & ' %s '\n",Name.version);
+fprintf(out,"      write(iout,*)'flags:'\n");
+}
+system("rm tmptmp");
+fclose(in1);
+in1=fopen("Makefile","r");
+while(fgets(buf,498,in1) != NULL)
+ {
+ if((p1=strchr(buf,'=')) != NULL && buf[0] != '#')
+  {
+  buf[strlen(buf)-1]='\0';
+  if(strlen(buf) > 49)
+   {
+   buf[47]='\0';
+   strcat(buf,"...");
+   }
+  else
+   {
+   while(buf[strlen(buf)-1]=='\\')
+    {
+    strcat(buf,"\\");
+    fprintf(out,"      write(iout,*)'%s'\n",buf);
+    if (fgets(buf,498,in1) != NULL)
+       buf[strlen(buf)-1]='\0';
+    if(strlen(buf) > 49)
+     {
+     buf[47]='\0';
+     strcat(buf,"...");
+     }
+    }
+   }
+  
+  fprintf(out,"      write(iout,*)'%s'\n",buf);
+  }
+ }
+fprintf(out,"      write(iout,*)'++++ End of compile info ++++'\n");
+fprintf(out,"      return\n");
+fprintf(out,"      end\n");
+fclose(out);
+fclose(in1);
+fclose(in);
+system("mv cinfo.f.new cinfo.f");
+}

diff --git a/source/unres/src_MIN/djacob.f b/source/unres/src_MIN/djacob.f
new file mode 100644 (file)
index 0000000..e3f46bc
--- /dev/null
+++ b/source/unres/src_MIN/djacob.f
@@ -0,0 +1,107 @@
+      SUBROUTINE DJACOB(N,NMAX,MAXJAC,E,A,C,AII)                        
+      IMPLICIT REAL*8 (A-H,O-Z) 
+C     THE JACOBI DIAGONALIZATION PROCEDURE
+      COMMON INP,IOUT,IPN                              
+      DIMENSION A(NMAX,N),C(NMAX,N),AII(150),AJJ(150)
+      SIN45 = .70710678                                                 
+      COS45 = .70710678                                                 
+      S45SQ = 0.50                                                      
+      C45SQ = 0.50                                                      
+C     UNIT EIGENVECTOR MATRIX                                           
+      DO 70 I = 1,N                                                     
+      DO 7 J = I,N                                                      
+      A(J,I)=A(I,J)                                                     
+      C(I,J) = 0.0                                                      
+    7 C(J,I) = 0.0                                                      
+   70  C(I,I) = 1.0                                                     
+C     DETERMINATION OF SEARCH ARGUMENT, TEST                            
+      AMAX = 0.0                                                        
+      DO 1 I = 1,N                                                      
+      DO 1 J = 1,I                                                      
+      TEMPA=DABS(A(I,J))                                                 
+      IF (AMAX-TEMPA) 2,1,1                                             
+    2 AMAX = TEMPA                                                      
+    1 CONTINUE                                                          
+      TEST = AMAX*E                                                     
+C     SEARCH FOR LARGEST OFF DIAGONAL ELEMENT                           
+      DO 72 IJAC=1,MAXJAC                                               
+      AIJMAX = 0.0                                                      
+      DO 3 I = 2,N                                                      
+      LIM = I-1                                                         
+      DO 3 J = 1,LIM                                                    
+      TAIJ=DABS(A(I,J))                                                  
+       IF (AIJMAX-TAIJ) 4,3,3                                           
+    4 AIJMAX = TAIJ                                                     
+      IPIV = I                                                          
+      JPIV = J                                                          
+    3 CONTINUE                                                          
+      IF(AIJMAX-TEST)300,300,5                                          
+C     PARAMETERS FOR ROTATION                                           
+    5 TAII = A(IPIV,IPIV)                                               
+      TAJJ = A(JPIV,JPIV)                                               
+      TAIJ = A(IPIV,JPIV)                                               
+      TMT = TAII-TAJJ                                                   
+      IF(DABS(TMT/TAIJ)-1.0D-12) 60,60,6                                 
+   60 IF(TAIJ) 10,10,11                                                 
+    6 ZAMMA=TAIJ/(2.0*TMT)                                              
+   90 IF(DABS(ZAMMA)-0.38268)8,8,9                                       
+    9 IF(ZAMMA)10,10,11                                                 
+   10 SINT = -SIN45                                                     
+      GO TO 12                                                          
+   11 SINT = SIN45                                                      
+   12 COST = COS45                                                      
+      SINSQ = S45SQ                                                     
+      COSSQ = C45SQ                                                     
+      GO TO 120                                                         
+    8 GAMSQ=ZAMMA*ZAMMA                                                 
+      SINT=2.0*ZAMMA/(1.0+GAMSQ)                                        
+      COST = (1.0-GAMSQ)/(1.0+GAMSQ)                                    
+      SINSQ=SINT*SINT                                                   
+      COSSQ=COST*COST                                                   
+C     ROTATION                                                          
+  120 DO 13 K = 1,N                                                     
+      TAIK = A(IPIV,K)                                                  
+      TAJK = A(JPIV,K)                                                  
+      A(IPIV,K) = TAIK*COST+TAJK*SINT                                   
+      A(JPIV,K) = TAJK*COST-TAIK*SINT                                   
+      TCIK = C(IPIV,K)                                                  
+      TCJK = C(JPIV,K)                                                  
+      C(IPIV,K) = TCIK*COST+TCJK*SINT                                   
+   13 C(JPIV,K) = TCJK*COST-TCIK*SINT                                   
+      A(IPIV,IPIV) = TAII*COSSQ+TAJJ*SINSQ+2.0*TAIJ*SINT*COST           
+      A(JPIV,JPIV) = TAII*SINSQ+TAJJ*COSSQ-2.0*TAIJ*SINT*COST           
+      A(IPIV,JPIV) = TAIJ*(COSSQ-SINSQ)-SINT*COST*TMT                   
+      A(JPIV,IPIV) = A(IPIV,JPIV)                                       
+      DO 30 K = 1,N                                                     
+      A(K,IPIV) = A(IPIV,K)                                             
+   30 A(K,JPIV) = A(JPIV,K)                                             
+   72 CONTINUE                                                          
+      WRITE (IOUT,1000) AIJMAX                                             
+ 1000 FORMAT (/1X,'NONCONVERGENT JACOBI. LARGEST OFF-DIAGONAL ELE',     
+     1 'MENT = ',1PE14.7)                                               
+C     ARRANGEMENT OF EIGENVALUES IN ASCENDING ORDER                     
+  300 DO 14 I=1,N                                                       
+   14 AJJ(I)=A(I,I)                                                     
+      LT=N+1                                                            
+      DO15 L=1,N                                                        
+      LT=LT-1                                                           
+      AIIMIN=1.0E+30                                                    
+      DO16 I=1,N                                                        
+      IF(AJJ(I)-AIIMIN)17,16,16                                         
+   17 AIIMIN=AJJ(I)                                                     
+      IT=I                                                              
+   16 CONTINUE                                                          
+      IN=L                                                              
+      AII(IN)=AIIMIN                                                    
+      AJJ(IT)=1.0E+30                                                   
+      DO15 K=1,N                                                        
+   15 A(IN,K)=C(IT,K)                                                   
+      DO 18 I=1,N                                                       
+      IF(A(I,1))19,22,22                                                
+   19 T=-1.0                                                            
+      GO TO 91                                                          
+   22 T=1.0                                                             
+   91 DO 18 J=1,N                                                       
+   18 C(J,I)=T*A(I,J)                                                   
+      RETURN                                                            
+      END

diff --git a/source/unres/src_MIN/gen_rand_conf.F b/source/unres/src_MIN/gen_rand_conf.F
new file mode 100644 (file)
index 0000000..6cc31ba
--- /dev/null
+++ b/source/unres/src_MIN/gen_rand_conf.F
@@ -0,0 +1,910 @@
+      subroutine gen_rand_conf(nstart,*)
+C Generate random conformation or chain cut and regrowth.
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.LOCAL'
+      include 'COMMON.VAR'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.MCM'
+      include 'COMMON.GEO'
+      include 'COMMON.CONTROL'
+      logical overlap,back,fail
+cd    print *,' CG Processor',me,' maxgen=',maxgen
+      maxsi=100
+cd    write (iout,*) 'Gen_Rand_conf: nstart=',nstart
+      if (nstart.lt.5) then
+        it1=itype(2)
+        phi(4)=gen_phi(4,itype(2),itype(3))
+c       write(iout,*)'phi(4)=',rad2deg*phi(4)
+        if (nstart.lt.3) theta(3)=gen_theta(itype(2),pi,phi(4))
+c       write(iout,*)'theta(3)=',rad2deg*theta(3) 
+        if (it1.ne.10) then
+          nsi=0
+          fail=.true.
+          do while (fail.and.nsi.le.maxsi)
+            call gen_side(it1,theta(3),alph(2),omeg(2),fail)
+            nsi=nsi+1
+          enddo
+          if (nsi.gt.maxsi) return1
+        endif ! it1.ne.10
+        call orig_frame
+        i=4
+        nstart=4
+      else
+        i=nstart
+        nstart=max0(i,4)
+      endif
+
+      maxnit=0
+
+      nit=0
+      niter=0
+      back=.false.
+      do while (i.le.nres .and. niter.lt.maxgen)
+        if (i.lt.nstart) then
+          if(iprint.gt.1) then
+          write (iout,'(/80(1h*)/2a/80(1h*))') 
+     &          'Generation procedure went down to ',
+     &          'chain beginning. Cannot continue...'
+          write (*,'(/80(1h*)/2a/80(1h*))') 
+     &          'Generation procedure went down to ',
+     &          'chain beginning. Cannot continue...'
+          endif
+          return1
+        endif
+       it1=itype(i-1)
+       it2=itype(i-2)
+       it=itype(i)
+c       print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
+c    &    ' nit=',nit,' niter=',niter,' maxgen=',maxgen
+       phi(i+1)=gen_phi(i+1,it1,it)
+       if (back) then
+          phi(i)=gen_phi(i+1,it2,it1)
+c         print *,'phi(',i,')=',phi(i)
+         theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
+         if (it2.ne.10) then
+            nsi=0
+            fail=.true.
+            do while (fail.and.nsi.le.maxsi)
+              call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail)
+              nsi=nsi+1
+            enddo
+            if (nsi.gt.maxsi) return1
+          endif
+         call locate_next_res(i-1)
+       endif
+       theta(i)=gen_theta(it1,phi(i),phi(i+1))
+        if (it1.ne.10) then 
+        nsi=0
+        fail=.true.
+        do while (fail.and.nsi.le.maxsi)
+          call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
+          nsi=nsi+1
+        enddo
+        if (nsi.gt.maxsi) return1
+        endif
+       call locate_next_res(i)
+       if (overlap(i-1)) then
+         if (nit.lt.maxnit) then
+           back=.true.
+           nit=nit+1
+          else
+           nit=0
+           if (i.gt.3) then
+             back=.true.
+             i=i-1
+            else
+             write (iout,'(a)') 
+     &  'Cannot generate non-overlaping conformation. Increase MAXNIT.'
+             write (*,'(a)') 
+     &  'Cannot generate non-overlaping conformation. Increase MAXNIT.'
+             return1
+           endif
+          endif
+        else
+         back=.false.
+         nit=0 
+         i=i+1
+        endif
+       niter=niter+1
+      enddo
+      if (niter.ge.maxgen) then
+       write (iout,'(a,2i5)') 
+     & 'Too many trials in conformation generation',niter,maxgen
+       write (*,'(a,2i5)') 
+     & 'Too many trials in conformation generation',niter,maxgen
+       return1
+      endif
+      do j=1,3
+       c(j,nres+1)=c(j,1)
+       c(j,nres+nres)=c(j,nres)
+      enddo
+      return
+      end
+c-------------------------------------------------------------------------
+      logical function overlap(i)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.FFIELD'
+      data redfac /0.5D0/
+      overlap=.false.
+      iti=itype(i)
+      if (iti.gt.ntyp) return
+C Check for SC-SC overlaps.
+cd    print *,'nnt=',nnt,' nct=',nct
+      do j=nnt,i-1
+        itj=itype(j)
+        if (j.lt.i-1 .or. ipot.ne.4) then
+          rcomp=sigmaii(iti,itj)
+        else 
+          rcomp=sigma(iti,itj)
+        endif
+cd      print *,'j=',j
+       if (dist(nres+i,nres+j).lt.redfac*rcomp) then
+          overlap=.true.
+c        print *,'overlap, SC-SC: i=',i,' j=',j,
+c     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
+c     &     rcomp
+         return
+        endif
+      enddo
+C Check for overlaps between the added peptide group and the preceding
+C SCs.
+      iteli=itel(i)
+      do j=1,3
+       c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
+      enddo
+      do j=nnt,i-2
+       itj=itype(j)
+cd      print *,'overlap, p-Sc: i=',i,' j=',j,
+cd   &         ' dist=',dist(nres+j,maxres2+1)
+       if (dist(nres+j,maxres2+1).lt.4.0D0*redfac) then
+         overlap=.true.
+         return
+        endif
+      enddo
+C Check for overlaps between the added side chain and the preceding peptide
+C groups.
+      do j=1,nnt-2
+       do k=1,3
+         c(k,maxres2+1)=0.5D0*(c(k,j)+c(k,j+1))
+        enddo
+cd      print *,'overlap, SC-p: i=',i,' j=',j,
+cd   &         ' dist=',dist(nres+i,maxres2+1)
+       if (dist(nres+i,maxres2+1).lt.4.0D0*redfac) then
+          overlap=.true.
+         return
+        endif
+      enddo
+C Check for p-p overlaps
+      do j=1,3
+       c(j,maxres2+2)=0.5D0*(c(j,i)+c(j,i+1))
+      enddo
+      do j=nnt,i-2
+        itelj=itel(j)
+       do k=1,3
+         c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1))
+        enddo
+cd      print *,'overlap, p-p: i=',i,' j=',j,
+cd   &         ' dist=',dist(maxres2+1,maxres2+2)
+        if(iteli.ne.0.and.itelj.ne.0)then
+        if (dist(maxres2+1,maxres2+2).lt.rpp(iteli,itelj)*redfac) then
+          overlap=.true.
+          return
+        endif
+        endif
+      enddo
+      return
+      end
+c--------------------------------------------------------------------------
+      double precision function gen_phi(i,it1,it2)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.GEO'
+      include 'COMMON.BOUNDS'
+c      gen_phi=ran_number(-pi,pi) 
+C 8/13/98 Generate phi using pre-defined boundaries
+      gen_phi=ran_number(phibound(1,i),phibound(2,i)) 
+      return
+      end
+c---------------------------------------------------------------------------
+      double precision function gen_theta(it,gama,gama1)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.LOCAL'
+      include 'COMMON.GEO'
+      double precision y(2),z(2)
+      double precision theta_max,theta_min
+c     print *,'gen_theta: it=',it
+      theta_min=0.05D0*pi
+      theta_max=0.95D0*pi
+      if (dabs(gama).gt.dwapi) then
+        y(1)=dcos(gama)
+        y(2)=dsin(gama)
+      else
+        y(1)=0.0D0
+        y(2)=0.0D0
+      endif
+      if (dabs(gama1).gt.dwapi) then
+        z(1)=dcos(gama1)
+        z(2)=dsin(gama1)
+      else
+       z(1)=0.0D0
+       z(2)=0.0D0
+      endif  
+      thet_pred_mean=a0thet(it)
+      do k=1,2
+        thet_pred_mean=thet_pred_mean+athet(k,it)*y(k)+bthet(k,it)*z(k)
+      enddo
+      sig=polthet(3,it)
+      do j=2,0,-1
+        sig=sig*thet_pred_mean+polthet(j,it)
+      enddo
+      sig=0.5D0/(sig*sig+sigc0(it))
+      ak=dexp(gthet(1,it)-
+     &0.5D0*((gthet(2,it)-thet_pred_mean)/gthet(3,it))**2)
+c     print '(i5,5(1pe14.4))',it,(gthet(j,it),j=1,3)
+c     print '(5(1pe14.4))',thet_pred_mean,theta0(it),sig,sig0(it),ak
+      theta_temp=binorm(thet_pred_mean,theta0(it),sig,sig0(it),ak) 
+      if (theta_temp.lt.theta_min) theta_temp=theta_min
+      if (theta_temp.gt.theta_max) theta_temp=theta_max
+      gen_theta=theta_temp
+c     print '(a)','Exiting GENTHETA.'
+      return
+      end
+c-------------------------------------------------------------------------
+      subroutine gen_side(it,the,al,om,fail)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.GEO'
+      include 'COMMON.LOCAL'
+      include 'COMMON.SETUP'
+      include 'COMMON.IOUNITS'
+      double precision MaxBoxLen /10.0D0/
+      double precision Ap_inv(3,3),a(3,3),z(3,maxlob),W1(maxlob),
+     & sumW(0:maxlob),y(2),cm(2),eig(2),box(2,2),work(100),detAp(maxlob)
+      double precision eig_limit /1.0D-8/
+      double precision Big /10.0D0/
+      double precision vec(3,3)
+      logical lprint,fail,lcheck
+      lcheck=.false.
+      lprint=.false.
+      fail=.false.
+      if (the.eq.0.0D0 .or. the.eq.pi) then
+#ifdef MPI
+        write (*,'(a,i4,a,i3,a,1pe14.5)') 
+     & 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the
+#else
+cd        write (iout,'(a,i3,a,1pe14.5)') 
+cd     &   'Error in GenSide: it=',it,' theta=',the
+#endif
+        fail=.true.
+        return
+      endif
+      tant=dtan(the-pipol)
+      nlobit=nlob(it)
+      if (lprint) then
+#ifdef MPI
+        print '(a,i4,a)','CG Processor:',me,' Enter Gen_Side.'
+        write (iout,'(a,i4,a)') 'Processor:',me,' Enter Gen_Side.'
+#endif
+        print *,'it=',it,' nlobit=',nlobit,' the=',the,' tant=',tant
+        write (iout,*) 'it=',it,' nlobit=',nlobit,' the=',the,
+     &     ' tant=',tant
+      endif
+      do i=1,nlobit
+       zz1=tant-censc(1,i,it)
+        do k=1,3
+          do l=1,3
+            a(k,l)=gaussc(k,l,i,it)
+          enddo
+        enddo
+        detApi=a(2,2)*a(3,3)-a(2,3)**2
+        Ap_inv(2,2)=a(3,3)/detApi
+        Ap_inv(2,3)=-a(2,3)/detApi
+        Ap_inv(3,2)=Ap_inv(2,3)
+        Ap_inv(3,3)=a(2,2)/detApi
+        if (lprint) then
+          write (*,'(/a,i2/)') 'Cluster #',i
+          write (*,'(3(1pe14.5),5x,1pe14.5)') 
+     &    ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
+          write (iout,'(/a,i2/)') 'Cluster #',i
+          write (iout,'(3(1pe14.5),5x,1pe14.5)') 
+     &    ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
+        endif
+        W1i=0.0D0
+        do k=2,3
+          do l=2,3
+            W1i=W1i+a(k,1)*a(l,1)*Ap_inv(k,l)
+          enddo
+        enddo
+        W1i=a(1,1)-W1i
+        W1(i)=dexp(bsc(i,it)-0.5D0*W1i*zz1*zz1)
+c        if (lprint) write(*,'(a,3(1pe15.5)/)')
+c     &          'detAp, W1, anormi',detApi,W1i,anormi
+       do k=2,3
+         zk=censc(k,i,it)
+         do l=2,3
+            zk=zk+zz1*Ap_inv(k,l)*a(l,1)
+          enddo
+         z(k,i)=zk
+        enddo
+        detAp(i)=dsqrt(detApi)
+      enddo
+
+      if (lprint) then
+        print *,'W1:',(w1(i),i=1,nlobit)
+        print *,'detAp:',(detAp(i),i=1,nlobit)
+        print *,'Z'
+        do i=1,nlobit
+          print '(i2,3f10.5)',i,(rad2deg*z(j,i),j=2,3)
+        enddo
+        write (iout,*) 'W1:',(w1(i),i=1,nlobit)
+        write (iout,*) 'detAp:',(detAp(i),i=1,nlobit)
+        write (iout,*) 'Z'
+        do i=1,nlobit
+          write (iout,'(i2,3f10.5)') i,(rad2deg*z(j,i),j=2,3)
+        enddo
+      endif
+      if (lcheck) then
+C Writing the distribution just to check the procedure
+      fac=0.0D0
+      dV=deg2rad**2*10.0D0
+      sum=0.0D0
+      sum1=0.0D0
+      do i=1,nlobit
+        fac=fac+W1(i)/detAp(i)
+      enddo 
+      fac=1.0D0/(2.0D0*fac*pi)
+cd    print *,it,'fac=',fac
+      do ial=90,180,2
+        y(1)=deg2rad*ial
+        do iom=-180,180,5
+          y(2)=deg2rad*iom
+          wart=0.0D0
+          do i=1,nlobit
+            do j=2,3
+              do k=2,3
+                a(j-1,k-1)=gaussc(j,k,i,it)
+              enddo
+            enddo
+            y2=y(2)
+
+            do iii=-1,1
+          
+              y(2)=y2+iii*dwapi
+
+              wykl=0.0D0
+              do j=1,2
+                do k=1,2 
+                  wykl=wykl+a(j,k)*(y(j)-z(j+1,i))*(y(k)-z(k+1,i))
+                enddo
+              enddo
+              wart=wart+W1(i)*dexp(-0.5D0*wykl)
+
+            enddo
+
+            y(2)=y2
+
+          enddo
+c         print *,'y',y(1),y(2),' fac=',fac
+          wart=fac*wart
+          write (20,'(2f10.3,1pd15.5)') y(1)*rad2deg,y(2)*rad2deg,wart
+          sum=sum+wart
+          sum1=sum1+1.0D0
+        enddo
+      enddo
+c     print *,'it=',it,' sum=',sum*dV,' sum1=',sum1*dV
+      return
+      endif
+
+C Calculate the CM of the system
+C
+      do i=1,nlobit
+        W1(i)=W1(i)/detAp(i)
+      enddo
+      sumW(0)=0.0D0
+      do i=1,nlobit
+       sumW(i)=sumW(i-1)+W1(i)
+      enddo
+      cm(1)=z(2,1)*W1(1)
+      cm(2)=z(3,1)*W1(1)
+      do j=2,nlobit
+        cm(1)=cm(1)+z(2,j)*W1(j) 
+        cm(2)=cm(2)+W1(j)*(z(3,1)+pinorm(z(3,j)-z(3,1)))
+      enddo
+      cm(1)=cm(1)/sumW(nlobit)
+      cm(2)=cm(2)/sumW(nlobit)
+      if (cm(1).gt.Big .or. cm(1).lt.-Big .or.
+     & cm(2).gt.Big .or. cm(2).lt.-Big) then
+cd        write (iout,'(a)') 
+cd     & 'Unexpected error in GenSide - CM coordinates too large.'
+cd        write (iout,'(i5,2(1pe14.5))') it,cm(1),cm(2)
+cd        write (*,'(a)') 
+cd     & 'Unexpected error in GenSide - CM coordinates too large.'
+cd        write (*,'(i5,2(1pe14.5))') it,cm(1),cm(2)
+        fail=.true. 
+        return
+      endif
+cd    print *,'CM:',cm(1),cm(2)
+C
+C Find the largest search distance from CM
+C
+      radmax=0.0D0
+      do i=1,nlobit
+       do j=2,3
+         do k=2,3
+           a(j-1,k-1)=gaussc(j,k,i,it) 
+          enddo
+       enddo
+#ifdef NAG
+        call f02faf('N','U',2,a,3,eig,work,100,ifail)
+#else
+        call djacob(2,3,10000,1.0d-10,a,vec,eig)
+#endif
+#ifdef MPI
+        if (lprint) then
+          print *,'*************** CG Processor',me
+          print *,'CM:',cm(1),cm(2)
+          write (iout,*) '*************** CG Processor',me
+          write (iout,*) 'CM:',cm(1),cm(2)
+          print '(A,8f10.5)','Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
+          write (iout,'(A,8f10.5)')
+     &        'Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
+        endif
+#endif
+        if (eig(1).lt.eig_limit) then
+          write(iout,'(a)')
+     &     'From Mult_Norm: Eigenvalues of A are too small.'
+          write(*,'(a)')
+     &     'From Mult_Norm: Eigenvalues of A are too small.'
+         fail=.true.
+          return
+        endif
+       radius=0.0D0
+cd      print *,'i=',i
+       do j=1,2
+         radius=radius+pinorm(z(j+1,i)-cm(j))**2
+        enddo
+       radius=dsqrt(radius)+3.0D0/dsqrt(eig(1))
+       if (radius.gt.radmax) radmax=radius
+      enddo
+      if (radmax.gt.pi) radmax=pi
+C
+C Determine the boundaries of the search rectangle.
+C
+      if (lprint) then
+        print '(a,4(1pe14.4))','W1: ',(W1(i),i=1,nlob(it) )
+        print '(a,4(1pe14.4))','radmax: ',radmax
+      endif
+      box(1,1)=dmax1(cm(1)-radmax,0.0D0)
+      box(2,1)=dmin1(cm(1)+radmax,pi)
+      box(1,2)=cm(2)-radmax
+      box(2,2)=cm(2)+radmax
+      if (lprint) then
+#ifdef MPI
+        print *,'CG Processor',me,' Array BOX:'
+#else
+        print *,'Array BOX:'
+#endif
+        print '(4(1pe14.4))',((box(k,j),k=1,2),j=1,2)
+        print '(a,4(1pe14.4))','sumW: ',(sumW(i),i=0,nlob(it) )
+#ifdef MPI
+        write (iout,*)'CG Processor',me,' Array BOX:'
+#else
+        write (iout,*)'Array BOX:'
+#endif
+        write(iout,'(4(1pe14.4))') ((box(k,j),k=1,2),j=1,2)
+        write(iout,'(a,4(1pe14.4))')'sumW: ',(sumW(i),i=0,nlob(it) )
+      endif
+      if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
+#ifdef MPI
+        write (iout,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
+        write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
+#else
+c        write (iout,'(a)') 'Bad sampling box.'
+#endif
+        fail=.true.
+        return
+      endif
+      which_lobe=ran_number(0.0D0,sumW(nlobit))
+c     print '(a,1pe14.4)','which_lobe=',which_lobe
+      do i=1,nlobit
+        if (sumW(i-1).le.which_lobe .and. sumW(i).ge.which_lobe) goto 1
+      enddo
+    1 ilob=i
+c     print *,'ilob=',ilob,' nlob=',nlob(it)
+      do i=2,3
+       cm(i-1)=z(i,ilob)
+       do j=2,3
+         a(i-1,j-1)=gaussc(i,j,ilob,it)
+        enddo
+      enddo
+cd    print '(a,i4,a)','CG Processor',me,' Calling MultNorm1.'
+      call mult_norm1(3,2,a,cm,box,y,fail)
+      if (fail) return
+      al=y(1)
+      om=pinorm(y(2))
+cd    print *,'al=',al,' om=',om
+cd    stop
+      return
+      end
+c---------------------------------------------------------------------------
+      double precision function ran_number(x1,x2)
+C Calculate a random real number from the range (x1,x2).
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      double precision x1,x2,fctor
+      data fctor /2147483647.0D0/
+#ifdef MPI
+      include "mpif.h"
+      include 'COMMON.SETUP'
+      ran_number=x1+(x2-x1)*prng_next(me)
+#else
+      call vrnd(ix,1)
+      ran_number=x1+(x2-x1)*ix/fctor
+#endif
+      return
+      end
+c--------------------------------------------------------------------------
+      integer function iran_num(n1,n2)
+C Calculate a random integer number from the range (n1,n2).
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      integer n1,n2,ix
+      real fctor /2147483647.0/
+#ifdef MPI
+      include "mpif.h"
+      include 'COMMON.SETUP'
+      ix=n1+(n2-n1+1)*prng_next(me)
+      if (ix.lt.n1) ix=n1
+      if (ix.gt.n2) ix=n2
+      iran_num=ix
+#else
+      call vrnd(ix,1)
+      ix=n1+(n2-n1+1)*(ix/fctor)
+      if (ix.gt.n2) ix=n2
+      iran_num=ix
+#endif
+      return
+      end
+c--------------------------------------------------------------------------
+      double precision function binorm(x1,x2,sigma1,sigma2,ak)
+      implicit real*8 (a-h,o-z)
+c     print '(a)','Enter BINORM.'
+      alowb=dmin1(x1-3.0D0*sigma1,x2-3.0D0*sigma2)
+      aupb=dmax1(x1+3.0D0*sigma1,x2+3.0D0*sigma2)
+      seg=sigma1/(sigma1+ak*sigma2)
+      alen=ran_number(0.0D0,1.0D0)
+      if (alen.lt.seg) then
+        binorm=anorm_distr(x1,sigma1,alowb,aupb)
+      else
+        binorm=anorm_distr(x2,sigma2,alowb,aupb)
+      endif
+c     print '(a)','Exiting BINORM.'
+      return
+      end
+c-----------------------------------------------------------------------
+c      double precision function anorm_distr(x,sigma,alowb,aupb)
+c      implicit real*8 (a-h,o-z)
+c     print '(a)','Enter ANORM_DISTR.'
+c   10 y=ran_number(alowb,aupb)
+c      expon=dexp(-0.5D0*((y-x)/sigma)**2)
+c      ran=ran_number(0.0D0,1.0D0)
+c      if (expon.lt.ran) goto 10
+c      anorm_distr=y
+c     print '(a)','Exiting ANORM_DISTR.'
+c      return
+c      end
+c-----------------------------------------------------------------------
+        double precision function anorm_distr(x,sigma,alowb,aupb)
+        implicit real*8 (a-h,o-z)
+c  to make a normally distributed deviate with zero mean and unit variance
+c
+        integer iset
+        real fac,gset,rsq,v1,v2,ran1
+        save iset,gset
+        data iset/0/
+        if(iset.eq.0) then
+1               v1=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
+                v2=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
+                rsq=v1**2+v2**2
+                if(rsq.ge.1.d0.or.rsq.eq.0.0d0) goto 1
+                fac=sqrt(-2.0d0*log(rsq)/rsq)
+                gset=v1*fac
+                gaussdev=v2*fac
+                iset=1
+        else
+                gaussdev=gset
+                iset=0
+        endif
+        anorm_distr=x+gaussdev*sigma
+      return
+      end
+c------------------------------------------------------------------------ 
+      subroutine mult_norm(lda,n,a,x,fail)
+C
+C Generate the vector X whose elements obey the multiple-normal distribution
+C from exp(-0.5*X'AX). LDA is the leading dimension of the moment matrix A,
+C n is the dimension of the problem. FAIL is set at .TRUE., if the smallest
+C eigenvalue of the matrix A is close to 0.
+C
+      implicit double precision (a-h,o-z)
+      double precision a(lda,n),x(n),eig(100),vec(3,3),work(100)
+      double precision eig_limit /1.0D-8/
+      logical fail
+      fail=.false.
+c     print '(a)','Enter MULT_NORM.'
+C
+C Find the smallest eigenvalue of the matrix A.
+C
+c     do i=1,n
+c       print '(8f10.5)',(a(i,j),j=1,n)
+c     enddo
+#ifdef NAG
+      call f02faf('V','U',2,a,lda,eig,work,100,ifail)
+#else
+      call djacob(2,lda,10000,1.0d-10,a,vec,eig)
+#endif
+c     print '(8f10.5)',(eig(i),i=1,n)
+C     print '(a)'
+c     do i=1,n
+c       print '(8f10.5)',(a(i,j),j=1,n)
+c     enddo
+      if (eig(1).lt.eig_limit) then
+        print *,'From Mult_Norm: Eigenvalues of A are too small.'
+        fail=.true.    
+       return
+      endif
+C 
+C Generate points following the normal distributions along the principal
+C axes of the moment matrix. Store in WORK.
+C
+      do i=1,n
+       sigma=1.0D0/dsqrt(eig(i))
+       alim=-3.0D0*sigma
+       work(i)=anorm_distr(0.0D0,sigma,-alim,alim)
+      enddo
+C
+C Transform the vector of normal variables back to the original basis.
+C
+      do i=1,n
+       xi=0.0D0
+       do j=1,n
+         xi=xi+a(i,j)*work(j)
+        enddo
+       x(i)=xi
+      enddo
+      return
+      end
+c------------------------------------------------------------------------ 
+      subroutine mult_norm1(lda,n,a,z,box,x,fail)
+C
+C Generate the vector X whose elements obey the multi-gaussian multi-dimensional
+C distribution from sum_{i=1}^m W(i)exp[-0.5*X'(i)A(i)X(i)]. LDA is the 
+C leading dimension of the moment matrix A, n is the dimension of the 
+C distribution, nlob is the number of lobes. FAIL is set at .TRUE., if the 
+C smallest eigenvalue of the matrix A is close to 0.
+C
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      double precision a(lda,n),z(n),x(n),box(n,n)
+      double precision etmp
+      include 'COMMON.IOUNITS'
+#ifdef MP
+      include 'COMMON.SETUP' 
+#endif
+      logical fail
+C 
+C Generate points following the normal distributions along the principal
+C axes of the moment matrix. Store in WORK.
+C
+cd    print *,'CG Processor',me,' entered MultNorm1.'
+cd    print '(2(1pe14.4),3x,1pe14.4)',((a(i,j),j=1,2),z(i),i=1,2)
+cd    do i=1,n
+cd      print *,i,box(1,i),box(2,i)
+cd    enddo
+      istep = 0
+   10 istep = istep + 1
+      if (istep.gt.10000) then
+c        write (iout,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
+c     & ' in MultNorm1.'
+c        write (*,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
+c     & ' in MultNorm1.'
+c        write (iout,*) 'box',box
+c        write (iout,*) 'a',a
+c        write (iout,*) 'z',z
+        fail=.true.
+        return
+      endif
+      do i=1,n
+       x(i)=ran_number(box(1,i),box(2,i))
+      enddo
+      ww=0.0D0
+      do i=1,n
+       xi=pinorm(x(i)-z(i))
+       ww=ww+0.5D0*a(i,i)*xi*xi
+       do j=i+1,n
+         ww=ww+a(i,j)*xi*pinorm(x(j)-z(j))
+        enddo
+      enddo
+      dec=ran_number(0.0D0,1.0D0)
+c      print *,(x(i),i=1,n),ww,dexp(-ww),dec
+crc   if (dec.gt.dexp(-ww)) goto 10
+      if(-ww.lt.100) then
+       etmp=dexp(-ww)
+      else
+       return  
+      endif
+      if (dec.gt.etmp) goto 10
+cd    print *,'CG Processor',me,' exitting MultNorm1.'
+      return
+      end
+c
+crc--------------------------------------
+      subroutine overlap_sc(scfail)
+c     Internal and cartesian coordinates must be consistent as input,
+c     and will be up-to-date on return.
+c     At the end of this procedure, scfail is true if there are
+c     overlapping residues left, or false otherwise (success)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.FFIELD'
+      include 'COMMON.VAR'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.IOUNITS'
+      logical had_overlaps,fail,scfail
+      integer ioverlap(maxres),ioverlap_last
+
+      had_overlaps=.false.
+      call overlap_sc_list(ioverlap,ioverlap_last)
+      if (ioverlap_last.gt.0) then
+        write (iout,*) '#OVERLAPing residues ',ioverlap_last
+        write (iout,'(20i4)') (ioverlap(k),k=1,ioverlap_last)
+        had_overlaps=.true.
+      endif
+
+      maxsi=1000
+      do k=1,1000
+        if (ioverlap_last.eq.0) exit
+
+        do ires=1,ioverlap_last 
+          i=ioverlap(ires)
+          iti=itype(i)
+          if (iti.ne.10) then
+            nsi=0
+            fail=.true.
+            do while (fail.and.nsi.le.maxsi)
+              call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
+              nsi=nsi+1
+            enddo
+            if(fail) goto 999
+          endif
+        enddo
+
+        call chainbuild
+        call overlap_sc_list(ioverlap,ioverlap_last)
+c        write (iout,*) 'Overlaping residues ',ioverlap_last,
+c     &           (ioverlap(j),j=1,ioverlap_last)
+      enddo
+
+      if (k.le.1000.and.ioverlap_last.eq.0) then
+        scfail=.false.
+        if (had_overlaps) then
+          write (iout,*) '#OVERLAPing all corrected after ',k,
+     &         ' random generation'
+        endif
+      else
+        scfail=.true.
+        write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
+        write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last)
+      endif
+
+      return
+
+ 999  continue
+      write (iout,'(a30,i5,a12,i4)') 
+     &               '#OVERLAP FAIL in gen_side after',maxsi,
+     &               'iter for RES',i
+      scfail=.true.
+      return
+      end
+
+      subroutine overlap_sc_list(ioverlap,ioverlap_last)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.GEO'
+      include 'COMMON.LOCAL'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.FFIELD'
+      include 'COMMON.VAR'
+      include 'COMMON.CALC'
+      logical fail
+      integer ioverlap(maxres),ioverlap_last
+      data redfac /0.5D0/
+
+      ioverlap_last=0
+C Check for SC-SC overlaps and mark residues
+c      print *,'>>overlap_sc nnt=',nnt,' nct=',nct
+      ind=0
+      do i=iatsc_s,iatsc_e
+        itypi=itype(i)
+        itypi1=itype(i+1)
+        xi=c(1,nres+i)
+        yi=c(2,nres+i)
+        zi=c(3,nres+i)
+        dxi=dc_norm(1,nres+i)
+        dyi=dc_norm(2,nres+i)
+        dzi=dc_norm(3,nres+i)
+        dsci_inv=dsc_inv(itypi)
+c
+       do iint=1,nint_gr(i)
+         do j=istart(i,iint),iend(i,iint)
+            ind=ind+1
+            itypj=itype(j)
+            dscj_inv=dsc_inv(itypj)
+            sig0ij=sigma(itypi,itypj)
+            chi1=chi(itypi,itypj)
+            chi2=chi(itypj,itypi)
+            chi12=chi1*chi2
+            chip1=chip(itypi)
+            chip2=chip(itypj)
+            chip12=chip1*chip2
+            alf1=alp(itypi)   
+            alf2=alp(itypj)   
+            alf12=0.5D0*(alf1+alf2)
+          if (j.gt.i+1) then
+           rcomp=sigmaii(itypi,itypj)
+          else 
+           rcomp=sigma(itypi,itypj)
+          endif
+c         print '(2(a3,2i3),a3,2f10.5)',
+c     &        ' i=',i,iti,' j=',j,itj,' d=',dist(nres+i,nres+j)
+c     &        ,rcomp
+            xj=c(1,nres+j)-xi
+            yj=c(2,nres+j)-yi
+            zj=c(3,nres+j)-zi
+            dxj=dc_norm(1,nres+j)
+            dyj=dc_norm(2,nres+j)
+            dzj=dc_norm(3,nres+j)
+            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+            rij=dsqrt(rrij)
+            call sc_angular
+            sigsq=1.0D0/sigsq
+            sig=sig0ij*dsqrt(sigsq)
+            rij_shift=1.0D0/rij-sig+sig0ij
+
+ct          if ( 1.0/rij .lt. redfac*rcomp .or. 
+ct     &       rij_shift.le.0.0D0 ) then
+            if ( rij_shift.le.0.0D0 ) then
+cd           write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)')
+cd     &     'overlap SC-SC: i=',i,' j=',j,
+cd     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
+cd     &     rcomp,1.0/rij,rij_shift
+          ioverlap_last=ioverlap_last+1
+          ioverlap(ioverlap_last)=i         
+          do k=1,ioverlap_last-1
+           if (ioverlap(k).eq.i) ioverlap_last=ioverlap_last-1
+          enddo
+          ioverlap_last=ioverlap_last+1
+          ioverlap(ioverlap_last)=j         
+          do k=1,ioverlap_last-1
+           if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1
+          enddo 
+         endif
+        enddo
+       enddo
+      enddo
+      return
+      end

diff --git a/source/unres/src_MIN/sc_move.F b/source/unres/src_MIN/sc_move.F
new file mode 100644 (file)
index 0000000..b6837fd
--- /dev/null
+++ b/source/unres/src_MIN/sc_move.F
@@ -0,0 +1,823 @@
+      subroutine sc_move(n_start,n_end,n_maxtry,e_drop,
+     +     n_fun,etot)
+c     Perform a quick search over side-chain arrangments (over
+c     residues n_start to n_end) for a given (frozen) CA trace
+c     Only side-chains are minimized (at most n_maxtry times each),
+c     not CA positions
+c     Stops if energy drops by e_drop, otherwise tries all residues
+c     in the given range
+c     If there is an energy drop, full minimization may be useful
+c     n_start, n_end CAN be modified by this routine, but only if
+c     out of bounds (n_start <= 1, n_end >= nres, n_start < n_end)
+c     NOTE: this move should never increase the energy
+crc      implicit none
+
+c     Includes
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      include 'COMMON.HEADER'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.FFIELD'
+
+c     External functions
+      integer iran_num
+      external iran_num
+
+c     Input arguments
+      integer n_start,n_end,n_maxtry
+      double precision e_drop
+
+c     Output arguments
+      integer n_fun
+      double precision etot
+
+c     Local variables
+      double precision energy(0:n_ene)
+      double precision cur_alph(2:nres-1),cur_omeg(2:nres-1)
+      double precision orig_e,cur_e
+      integer n,n_steps,n_first,n_cur,n_tot,i
+      double precision orig_w(n_ene)
+      double precision wtime
+
+
+c     Set non side-chain weights to zero (minimization is faster)
+c     NOTE: e(2) does not actually depend on the side-chain, only CA
+      orig_w(2)=wscp
+      orig_w(3)=welec
+      orig_w(4)=wcorr
+      orig_w(5)=wcorr5
+      orig_w(6)=wcorr6
+      orig_w(7)=wel_loc
+      orig_w(8)=wturn3
+      orig_w(9)=wturn4
+      orig_w(10)=wturn6
+      orig_w(11)=wang
+      orig_w(13)=wtor
+      orig_w(14)=wtor_d
+      orig_w(15)=wvdwpp
+
+      wscp=0.D0
+      welec=0.D0
+      wcorr=0.D0
+      wcorr5=0.D0
+      wcorr6=0.D0
+      wel_loc=0.D0
+      wturn3=0.D0
+      wturn4=0.D0
+      wturn6=0.D0
+      wang=0.D0
+      wtor=0.D0
+      wtor_d=0.D0
+      wvdwpp=0.D0
+
+c     Make sure n_start, n_end are within proper range
+      if (n_start.lt.2) n_start=2
+      if (n_end.gt.nres-1) n_end=nres-1
+crc      if (n_start.lt.n_end) then
+      if (n_start.gt.n_end) then
+        n_start=2
+        n_end=nres-1
+      endif
+
+c     Save the initial values of energy and coordinates
+cd      call chainbuild
+cd      call etotal(energy)
+cd      write (iout,*) 'start sc ene',energy(0)
+cd      call enerprint(energy(0))
+crc      etot=energy(0)
+       n_fun=0
+crc      orig_e=etot
+crc      cur_e=orig_e
+crc      do i=2,nres-1
+crc        cur_alph(i)=alph(i)
+crc        cur_omeg(i)=omeg(i)
+crc      enddo
+
+ct      wtime=MPI_WTIME()
+c     Try (one by one) all specified residues, starting from a
+c     random position in sequence
+c     Stop early if the energy has decreased by at least e_drop
+      n_tot=n_end-n_start+1
+      n_first=iran_num(0,n_tot-1)
+      n_steps=0
+      n=0
+crc      do while (n.lt.n_tot .and. orig_e-etot.lt.e_drop)
+      do while (n.lt.n_tot)
+        n_cur=n_start+mod(n_first+n,n_tot)
+        call single_sc_move(n_cur,n_maxtry,e_drop,
+     +       n_steps,n_fun,etot)
+c     If a lower energy was found, update the current structure...
+crc        if (etot.lt.cur_e) then
+crc          cur_e=etot
+crc          do i=2,nres-1
+crc            cur_alph(i)=alph(i)
+crc            cur_omeg(i)=omeg(i)
+crc          enddo
+crc        else
+c     ...else revert to the previous one
+crc          etot=cur_e
+crc          do i=2,nres-1
+crc            alph(i)=cur_alph(i)
+crc            omeg(i)=cur_omeg(i)
+crc          enddo
+crc        endif
+        n=n+1
+cd
+cd      call chainbuild
+cd      call etotal(energy)
+cd      print *,'running',n,energy(0)
+      enddo
+
+cd      call chainbuild
+cd      call etotal(energy)
+cd      write (iout,*) 'end   sc ene',energy(0)
+
+c     Put the original weights back to calculate the full energy
+      wscp=orig_w(2)
+      welec=orig_w(3)
+      wcorr=orig_w(4)
+      wcorr5=orig_w(5)
+      wcorr6=orig_w(6)
+      wel_loc=orig_w(7)
+      wturn3=orig_w(8)
+      wturn4=orig_w(9)
+      wturn6=orig_w(10)
+      wang=orig_w(11)
+      wtor=orig_w(13)
+      wtor_d=orig_w(14)
+      wvdwpp=orig_w(15)
+
+crc      n_fun=n_fun+1
+ct      write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime
+      return
+      end
+
+c-------------------------------------------------------------
+
+      subroutine single_sc_move(res_pick,n_maxtry,e_drop,
+     +     n_steps,n_fun,e_sc)
+c     Perturb one side-chain (res_pick) and minimize the
+c     neighbouring region, keeping all CA's and non-neighbouring
+c     side-chains fixed
+c     Try until e_drop energy improvement is achieved, or n_maxtry
+c     attempts have been made
+c     At the start, e_sc should contain the side-chain-only energy(0)
+c     nsteps and nfun for this move are ADDED to n_steps and n_fun
+crc      implicit none
+
+c     Includes
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.VAR'
+      include 'COMMON.INTERACT'
+      include 'COMMON.CHAIN'
+      include 'COMMON.MINIM'
+      include 'COMMON.FFIELD'
+      include 'COMMON.IOUNITS'
+
+c     External functions
+      double precision dist
+      external dist
+
+c     Input arguments
+      integer res_pick,n_maxtry
+      double precision e_drop
+
+c     Input/Output arguments
+      integer n_steps,n_fun
+      double precision e_sc
+
+c     Local variables
+      logical fail
+      integer i,j
+      integer nres_moved
+      integer iretcode,loc_nfun,orig_maxfun,n_try
+      double precision sc_dist,sc_dist_cutoff
+      double precision energy(0:n_ene),orig_e,cur_e
+      double precision evdw,escloc
+      double precision cur_alph(2:nres-1),cur_omeg(2:nres-1)
+      double precision var(maxvar)
+
+      double precision orig_theta(1:nres),orig_phi(1:nres),
+     +     orig_alph(1:nres),orig_omeg(1:nres)
+
+
+c     Define what is meant by "neighbouring side-chain"
+      sc_dist_cutoff=5.0D0
+
+c     Don't do glycine or ends
+      i=itype(res_pick)
+      if (i.eq.10 .or. i.eq.21) return
+
+c     Freeze everything (later will relax only selected side-chains)
+      mask_r=.true.
+      do i=1,nres
+        mask_phi(i)=0
+        mask_theta(i)=0
+        mask_side(i)=0
+      enddo
+
+c     Find the neighbours of the side-chain to move
+c     and save initial variables
+crc      orig_e=e_sc
+crc      cur_e=orig_e
+      nres_moved=0
+      do i=2,nres-1
+c     Don't do glycine (itype(j)==10)
+        if (itype(i).ne.10) then
+          sc_dist=dist(nres+i,nres+res_pick)
+        else
+          sc_dist=sc_dist_cutoff
+        endif
+        if (sc_dist.lt.sc_dist_cutoff) then
+          nres_moved=nres_moved+1
+          mask_side(i)=1
+          cur_alph(i)=alph(i)
+          cur_omeg(i)=omeg(i)
+        endif
+      enddo
+
+      call chainbuild
+      call egb1(evdw)
+      call esc(escloc)
+      e_sc=wsc*evdw+wscloc*escloc
+cd      call etotal(energy)
+cd      print *,'new       ',(energy(k),k=0,n_ene)
+      orig_e=e_sc
+      cur_e=orig_e
+
+      n_try=0
+      do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop)
+c     Move the selected residue (don't worry if it fails)
+        call gen_side(itype(res_pick),theta(res_pick+1),
+     +       alph(res_pick),omeg(res_pick),fail)
+
+c     Minimize the side-chains starting from the new arrangement
+        call geom_to_var(nvar,var)
+        orig_maxfun=maxfun
+        maxfun=7
+
+crc        do i=1,nres
+crc          orig_theta(i)=theta(i)
+crc          orig_phi(i)=phi(i)
+crc          orig_alph(i)=alph(i)
+crc          orig_omeg(i)=omeg(i)
+crc        enddo
+
+        call minimize_sc1(e_sc,var,iretcode,loc_nfun)
+        
+cv        write(*,'(2i3,2f12.5,2i3)') 
+cv     &       res_pick,nres_moved,orig_e,e_sc-cur_e,
+cv     &       iretcode,loc_nfun
+
+c$$$        if (iretcode.eq.8) then
+c$$$          write(iout,*)'Coordinates just after code 8'
+c$$$          call chainbuild
+c$$$          call all_varout
+c$$$          call flush(iout)
+c$$$          do i=1,nres
+c$$$            theta(i)=orig_theta(i)
+c$$$            phi(i)=orig_phi(i)
+c$$$            alph(i)=orig_alph(i)
+c$$$            omeg(i)=orig_omeg(i)
+c$$$          enddo
+c$$$          write(iout,*)'Coordinates just before code 8'
+c$$$          call chainbuild
+c$$$          call all_varout
+c$$$          call flush(iout)
+c$$$        endif
+
+        n_fun=n_fun+loc_nfun
+        maxfun=orig_maxfun
+        call var_to_geom(nvar,var)
+
+c     If a lower energy was found, update the current structure...
+        if (e_sc.lt.cur_e) then
+cv              call chainbuild
+cv              call etotal(energy)
+cd              call egb1(evdw)
+cd              call esc(escloc)
+cd              e_sc1=wsc*evdw+wscloc*escloc
+cd              print *,'     new',e_sc1,energy(0)
+cv              print *,'new       ',energy(0)
+cd              call enerprint(energy(0))
+          cur_e=e_sc
+          do i=2,nres-1
+            if (mask_side(i).eq.1) then
+              cur_alph(i)=alph(i)
+              cur_omeg(i)=omeg(i)
+            endif
+          enddo
+        else
+c     ...else revert to the previous one
+          e_sc=cur_e
+          do i=2,nres-1
+            if (mask_side(i).eq.1) then
+              alph(i)=cur_alph(i)
+              omeg(i)=cur_omeg(i)
+            endif
+          enddo
+        endif
+        n_try=n_try+1
+
+      enddo
+      n_steps=n_steps+n_try
+
+c     Reset the minimization mask_r to false
+      mask_r=.false.
+
+      return
+      end
+
+c-------------------------------------------------------------
+
+      subroutine sc_minimize(etot,iretcode,nfun)
+c     Minimizes side-chains only, leaving backbone frozen
+crc      implicit none
+
+c     Includes
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.VAR'
+      include 'COMMON.CHAIN'
+      include 'COMMON.FFIELD'
+
+c     Output arguments
+      double precision etot
+      integer iretcode,nfun
+
+c     Local variables
+      integer i
+      double precision orig_w(n_ene),energy(0:n_ene)
+      double precision var(maxvar)
+
+
+c     Set non side-chain weights to zero (minimization is faster)
+c     NOTE: e(2) does not actually depend on the side-chain, only CA
+      orig_w(2)=wscp
+      orig_w(3)=welec
+      orig_w(4)=wcorr
+      orig_w(5)=wcorr5
+      orig_w(6)=wcorr6
+      orig_w(7)=wel_loc
+      orig_w(8)=wturn3
+      orig_w(9)=wturn4
+      orig_w(10)=wturn6
+      orig_w(11)=wang
+      orig_w(13)=wtor
+      orig_w(14)=wtor_d
+
+      wscp=0.D0
+      welec=0.D0
+      wcorr=0.D0
+      wcorr5=0.D0
+      wcorr6=0.D0
+      wel_loc=0.D0
+      wturn3=0.D0
+      wturn4=0.D0
+      wturn6=0.D0
+      wang=0.D0
+      wtor=0.D0
+      wtor_d=0.D0
+
+c     Prepare to freeze backbone
+      do i=1,nres
+        mask_phi(i)=0
+        mask_theta(i)=0
+        mask_side(i)=1
+      enddo
+
+c     Minimize the side-chains
+      mask_r=.true.
+      call geom_to_var(nvar,var)
+      call minimize(etot,var,iretcode,nfun)
+      call var_to_geom(nvar,var)
+      mask_r=.false.
+
+c     Put the original weights back and calculate the full energy
+      wscp=orig_w(2)
+      welec=orig_w(3)
+      wcorr=orig_w(4)
+      wcorr5=orig_w(5)
+      wcorr6=orig_w(6)
+      wel_loc=orig_w(7)
+      wturn3=orig_w(8)
+      wturn4=orig_w(9)
+      wturn6=orig_w(10)
+      wang=orig_w(11)
+      wtor=orig_w(13)
+      wtor_d=orig_w(14)
+
+      call chainbuild
+      call etotal(energy)
+      etot=energy(0)
+
+      return
+      end
+
+c-------------------------------------------------------------
+      subroutine minimize_sc1(etot,x,iretcode,nfun)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
+      include 'COMMON.IOUNITS'
+      include 'COMMON.VAR'
+      include 'COMMON.GEO'
+      include 'COMMON.MINIM'
+      common /srutu/ icall
+      dimension iv(liv)                                               
+      double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
+      double precision energia(0:n_ene)
+      external func,gradient,fdum
+      external func_restr1,grad_restr1
+      logical not_done,change,reduce 
+      common /przechowalnia/ v
+
+      call deflt(2,iv,liv,lv,v)                                         
+* 12 means fresh start, dont call deflt                                 
+      iv(1)=12                                                          
+* max num of fun calls                                                  
+      if (maxfun.eq.0) maxfun=500
+      iv(17)=maxfun
+* max num of iterations                                                 
+      if (maxmin.eq.0) maxmin=1000
+      iv(18)=maxmin
+* controls output                                                       
+      iv(19)=2                                                          
+* selects output unit                                                   
+c     iv(21)=iout                                                       
+      iv(21)=0
+* 1 means to print out result                                           
+      iv(22)=0                                                          
+* 1 means to print out summary stats                                    
+      iv(23)=0                                                          
+* 1 means to print initial x and d                                      
+      iv(24)=0                                                          
+* min val for v(radfac) default is 0.1                                  
+      v(24)=0.1D0                                                       
+* max val for v(radfac) default is 4.0                                  
+      v(25)=2.0D0                                                       
+c     v(25)=4.0D0                                                       
+* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
+* the sumsl default is 0.1                                              
+      v(26)=0.1D0
+* false conv if (act fnctn decrease) .lt. v(34)                         
+* the sumsl default is 100*machep                                       
+      v(34)=v(34)/100.0D0                                               
+* absolute convergence                                                  
+      if (tolf.eq.0.0D0) tolf=1.0D-4
+      v(31)=tolf
+* relative convergence                                                  
+      if (rtolf.eq.0.0D0) rtolf=1.0D-4
+      v(32)=rtolf
+* controls initial step size                                            
+       v(35)=1.0D-1                                                    
+* large vals of d correspond to small components of step                
+      do i=1,nphi
+        d(i)=1.0D-1
+      enddo
+      do i=nphi+1,nvar
+        d(i)=1.0D-1
+      enddo
+      IF (mask_r) THEN
+       call x2xx(x,xx,nvar_restr)
+       call sumsl(nvar_restr,d,xx,func_restr1,grad_restr1,
+     &                    iv,liv,lv,v,idum,rdum,fdum)      
+       call xx2x(x,xx)
+      ELSE
+       call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)      
+      ENDIF
+      etot=v(10)                                                      
+      iretcode=iv(1)
+      nfun=iv(6)
+
+      return  
+      end  
+************************************************************************
+      subroutine func_restr1(n,x,nf,f,uiparm,urparm,ufparm)  
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.DERIV'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.GEO'
+      include 'COMMON.FFIELD'
+      include 'COMMON.INTERACT'
+      include 'COMMON.TIME1'
+      common /chuju/ jjj
+      double precision energia(0:n_ene),evdw,escloc
+      integer jjj
+      double precision ufparm,e1,e2
+      external ufparm                                                   
+      integer uiparm(1)                                                 
+      real*8 urparm(1)                                                    
+      dimension x(maxvar)
+      nfl=nf
+      icg=mod(nf,2)+1
+
+#ifdef OSF
+c     Intercept NaNs in the coordinates, before calling etotal
+      x_sum=0.D0
+      do i=1,n
+        x_sum=x_sum+x(i)
+      enddo
+      FOUND_NAN=.false.
+      if (x_sum.ne.x_sum) then
+        write(iout,*)"   *** func_restr1 : Found NaN in coordinates"
+        f=1.0D+73
+        FOUND_NAN=.true.
+        return
+      endif
+#endif
+
+      call var_to_geom_restr(n,x)
+      call zerograd
+      call chainbuild
+cd    write (iout,*) 'ETOTAL called from FUNC'
+      call egb1(evdw)
+      call esc(escloc)
+      f=wsc*evdw+wscloc*escloc
+cd      call etotal(energia(0))
+cd      f=wsc*energia(1)+wscloc*energia(12)
+cd      print *,f,evdw,escloc,energia(0)
+C
+C Sum up the components of the Cartesian gradient.
+C
+      do i=1,nct
+        do j=1,3
+          gradx(j,i,icg)=wsc*gvdwx(j,i)
+        enddo
+      enddo
+
+      return                                                            
+      end                                                               
+c-------------------------------------------------------
+      subroutine grad_restr1(n,x,nf,g,uiparm,urparm,ufparm)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.VAR'
+      include 'COMMON.INTERACT'
+      include 'COMMON.FFIELD'
+      include 'COMMON.IOUNITS'
+      external ufparm
+      integer uiparm(1)
+      double precision urparm(1)
+      dimension x(maxvar),g(maxvar)
+
+      icg=mod(nf,2)+1
+      if (nf-nfl+1) 20,30,40
+   20 call func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
+c     write (iout,*) 'grad 20'
+      if (nf.eq.0) return
+      goto 40
+   30 call var_to_geom_restr(n,x)
+      call chainbuild 
+C
+C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
+C
+   40 call cartder
+C
+C Convert the Cartesian gradient into internal-coordinate gradient.
+C
+
+      ig=0
+      ind=nres-2                                                                    
+      do i=2,nres-2                
+       IF (mask_phi(i+2).eq.1) THEN                                             
+        gphii=0.0D0                                                             
+        do j=i+1,nres-1                                                         
+          ind=ind+1                                 
+          do k=1,3                                                              
+            gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)                            
+            gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)                           
+          enddo                                                                 
+        enddo                                                                   
+        ig=ig+1
+        g(ig)=gphii
+       ELSE
+        ind=ind+nres-1-i
+       ENDIF
+      enddo                                        
+
+
+      ind=0
+      do i=1,nres-2
+       IF (mask_theta(i+2).eq.1) THEN
+        ig=ig+1
+       gthetai=0.0D0
+       do j=i+1,nres-1
+          ind=ind+1
+         do k=1,3
+            gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
+            gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
+          enddo
+        enddo
+        g(ig)=gthetai
+       ELSE
+        ind=ind+nres-1-i
+       ENDIF
+      enddo
+
+      do i=2,nres-1
+       if (itype(i).ne.10) then
+         IF (mask_side(i).eq.1) THEN
+          ig=ig+1
+          galphai=0.0D0
+         do k=1,3
+           galphai=galphai+dxds(k,i)*gradx(k,i,icg)
+          enddo
+          g(ig)=galphai
+         ENDIF
+        endif
+      enddo
+
+      
+      do i=2,nres-1
+        if (itype(i).ne.10) then
+         IF (mask_side(i).eq.1) THEN
+          ig=ig+1
+         gomegai=0.0D0
+         do k=1,3
+           gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
+          enddo
+         g(ig)=gomegai
+         ENDIF
+        endif
+      enddo
+
+C
+C Add the components corresponding to local energy terms.
+C
+
+      ig=0
+      igall=0
+      do i=4,nres
+        igall=igall+1
+        if (mask_phi(i).eq.1) then
+          ig=ig+1
+          g(ig)=g(ig)+gloc(igall,icg)
+        endif
+      enddo
+
+      do i=3,nres
+        igall=igall+1
+        if (mask_theta(i).eq.1) then
+          ig=ig+1
+          g(ig)=g(ig)+gloc(igall,icg)
+        endif
+      enddo
+     
+      do ij=1,2
+      do i=2,nres-1
+        if (itype(i).ne.10) then
+          igall=igall+1
+          if (mask_side(i).eq.1) then
+            ig=ig+1
+            g(ig)=g(ig)+gloc(igall,icg)
+          endif
+        endif
+      enddo
+      enddo
+
+cd      do i=1,ig
+cd        write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
+cd      enddo
+      return
+      end
+C-----------------------------------------------------------------------------
+      subroutine egb1(evdw)
+C
+C This subroutine calculates the interaction energy of nonbonded side chains
+C assuming the Gay-Berne potential of interaction.
+C
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      include 'COMMON.LOCAL'
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.NAMES'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CALC'
+      include 'COMMON.CONTROL'
+      logical lprn
+      evdw=0.0D0
+c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
+      evdw=0.0D0
+      lprn=.false.
+c     if (icall.eq.0) lprn=.true.
+      ind=0
+      do i=iatsc_s,iatsc_e
+
+
+        itypi=itype(i)
+        itypi1=itype(i+1)
+        xi=c(1,nres+i)
+        yi=c(2,nres+i)
+        zi=c(3,nres+i)
+        dxi=dc_norm(1,nres+i)
+        dyi=dc_norm(2,nres+i)
+        dzi=dc_norm(3,nres+i)
+        dsci_inv=dsc_inv(itypi)
+C
+C Calculate SC interaction energy.
+C
+        do iint=1,nint_gr(i)
+          do j=istart(i,iint),iend(i,iint)
+          IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
+            ind=ind+1
+            itypj=itype(j)
+            dscj_inv=dsc_inv(itypj)
+            sig0ij=sigma(itypi,itypj)
+            chi1=chi(itypi,itypj)
+            chi2=chi(itypj,itypi)
+            chi12=chi1*chi2
+            chip1=chip(itypi)
+            chip2=chip(itypj)
+            chip12=chip1*chip2
+            alf1=alp(itypi)
+            alf2=alp(itypj)
+            alf12=0.5D0*(alf1+alf2)
+C For diagnostics only!!!
+c           chi1=0.0D0
+c           chi2=0.0D0
+c           chi12=0.0D0
+c           chip1=0.0D0
+c           chip2=0.0D0
+c           chip12=0.0D0
+c           alf1=0.0D0
+c           alf2=0.0D0
+c           alf12=0.0D0
+            xj=c(1,nres+j)-xi
+            yj=c(2,nres+j)-yi
+            zj=c(3,nres+j)-zi
+            dxj=dc_norm(1,nres+j)
+            dyj=dc_norm(2,nres+j)
+            dzj=dc_norm(3,nres+j)
+            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+            rij=dsqrt(rrij)
+C Calculate angle-dependent terms of energy and contributions to their
+C derivatives.
+            call sc_angular
+            sigsq=1.0D0/sigsq
+            sig=sig0ij*dsqrt(sigsq)
+            rij_shift=1.0D0/rij-sig+sig0ij
+C I hate to put IF's in the loops, but here don't have another choice!!!!
+            if (rij_shift.le.0.0D0) then
+              evdw=1.0D20
+cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+cd     &        restyp(itypi),i,restyp(itypj),j,
+cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
+              return
+            endif
+            sigder=-sig*sigsq
+c---------------------------------------------------------------
+            rij_shift=1.0D0/rij_shift 
+            fac=rij_shift**expon
+            e1=fac*fac*aa(itypi,itypj)
+            e2=fac*bb(itypi,itypj)
+            evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+            eps2der=evdwij*eps3rt
+            eps3der=evdwij*eps2rt
+            evdwij=evdwij*eps2rt*eps3rt
+            evdw=evdw+evdwij
+            if (lprn) then
+            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+cd            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+cd     &        restyp(itypi),i,restyp(itypj),j,
+cd     &        epsi,sigm,chi1,chi2,chip1,chip2,
+cd     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+cd     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
+cd     &        evdwij
+            endif
+
+            if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') 
+     &                        'evdw',i,j,evdwij
+
+C Calculate gradient components.
+            e1=e1*eps1*eps2rt**2*eps3rt**2
+            fac=-expon*(e1+evdwij)*rij_shift
+            sigder=fac*sigder
+            fac=rij*fac
+C Calculate the radial part of the gradient
+            gg(1)=xj*fac
+            gg(2)=yj*fac
+            gg(3)=zj*fac
+C Calculate angular part of the gradient.
+            call sc_grad
+          ENDIF
+          enddo      ! j
+        enddo        ! iint
+      enddo          ! i
+      end
+C-----------------------------------------------------------------------------