Naprawienie wham'a oraz wprowadzenie SCcorr do CSA
authorAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Sun, 4 Jan 2015 10:54:16 +0000 (11:54 +0100)
committerAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Sun, 4 Jan 2015 10:54:16 +0000 (11:54 +0100)
19 files changed:
examples/unres/CSA/4P/1l2y_csa.inp
examples/unres/CSA/4P/start.mat
source/unres/src_CSA/CMakeLists.txt
source/unres/src_CSA/COMMON.CHAIN
source/unres/src_CSA/COMMON.LOCAL
source/unres/src_CSA/COMMON.SCCOR
source/unres/src_CSA/COMMON.VAR
source/unres/src_CSA/DIMENSIONS
source/unres/src_CSA/checkder_p.F
source/unres/src_CSA/energy_p_new_barrier.F
source/unres/src_CSA/gradient_p.F
source/unres/src_CSA/initialize_p.F
source/unres/src_CSA/int_to_cart.f
source/unres/src_CSA/intcartderiv.F
source/unres/src_CSA/parmread.F
source/unres/src_CSA/readpdb.F
source/unres/src_MD-M/int_to_cart.f
source/wham/src-M/COMMON.ALLPARM
source/wham/src-M/include_unres/COMMON.LOCAL

index 5c6fc33..f66f1ae 100644 (file)
@@ -8,7 +8,7 @@ ICMAX=1 IUCUT=2 NBANKM=200
 WLONG=1.00000 WSCP=2.73684 WELEC=0.06833 WANG=4.15526 WSCLOC=0.16761           &
 WTOR=2.99546 WTORD=2.89720 WCORRH=1.98989 WCORR5=0.00000 WCORR6=0.00000        &
 WEL_LOC=1.60072 WTURN3=2.36351 WTURN4=1.34051 WTURN6=0.00000                   &
-WSCCOR=0.0                                                                     &
+WSCCOR=0.1                                                                     &
 CUTOFF=7.00000 WCORR4=0.00000                                                  
 22
  D   ASN LEU TYR ILE GLN TRP LEU LYS ASP GLY GLY PRO SER SER GLY ARG PRO PRO PRO 
index bb5b966..e4e4073 100755 (executable)
@@ -1,9 +1,9 @@
 #PBS -N global_1l2y_csa
 #PBS -q special
-#PBS -l nodes=8:ppn=4
+#PBS -l nodes=2:ppn=4
 #PBS -l walltime=168:00:00
 
-setenv MPIRUN "/users/software/mpich2-1.4.1p1_gnu/bin/mpirun "
+setenv MPIRUN "/users/software/mpich2-1.4.1p1_intel/bin/mpirun "
 setenv NPROCS `cat $PBS_NODEFILE | wc -l`
 
 cd $PBS_O_WORKDIR
@@ -13,8 +13,8 @@ setenv POT GB
 setenv PREFIX 1l2y_csa
 setenv OUT1FILE NO
 #-----------------------------------------------------------------------------
-setenv UNRES_BIN /users2/czarek/UNRES/GIT/unres/bin/unres/CSA/unres_csa_gfort_MPICH_4P.exe
-setenv LD_LIBRARY_PATH /users/software/mpich2-1.4.1p1_gnu/lib:$LD_LIBRARY_PATH
+setenv UNRES_BIN /users2/adasko/unres/build/bin/unresCSA_ifort_4P.exe
+setenv LD_LIBRARY_PATH /users/software/mpich2-1.4.1p1_intel/lib:$LD_LIBRARY_PATH
 #-----------------------------------------------------------------------------
 setenv DD /users2/czarek/UNRES/GIT/unres/PARAM 
 setenv BONDPAR $DD/bond.parm
@@ -25,7 +25,8 @@ setenv TORDPAR $DD/torsion_double_631Gdp.parm
 setenv ELEPAR $DD/electr_631Gdp.parm
 setenv SIDEPAR $DD/sc_GB_opt.4P5_iter33_3r
 setenv FOURIER $DD/fourier_opt.parm.1igd_hc_iter3_3
-setenv SCCORPAR $DD/rotcorr_AM1.parm
+#setenv SCCORPAR $DD/rotcorr_AM1.parm
+setenv SCCORPAR $DD/sccor_am1_pawel.dat
 setenv SCPPAR $DD/scp.parm
 setenv PATTERN $DD/patterns.cart
 setenv PRINT_PARM NO
index 1710157..fdea19f 100644 (file)
@@ -148,7 +148,7 @@ elseif(UNRES_CSA_FF STREQUAL "CASP5")
 elseif(UNRES_CSA_FF STREQUAL "3P")
   set(CPPFLAGS "UNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC" )
 elseif(UNRES_CSA_FF STREQUAL "4P")
-  set(CPPFLAGS "UNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC" )
+  set(CPPFLAGS "UNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DSCCORPDB" )
 endif(UNRES_CSA_FF STREQUAL "CASP3")
 
 #=========================================
index f7a8a1d..6e19f8d 100644 (file)
@@ -1,9 +1,10 @@
       integer nres,nsup,nstart_sup,nz_start,nz_end,iz_sc,
      &  nres0,nstart_seq
       double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm,t,r,
-     & prod,rt,dc_work,cref,crefjlee
+     & prod,rt,dc_work,cref,crefjlee,dc_norm2
       common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
      & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
+     & dc_norm2(3,0:maxres2),
      & dc_work(MAXRES6),nres,nres0
       common /rotmat/ t(3,3,maxres),r(3,3,maxres),prod(3,3,maxres),
      &                rt(3,3,maxres) 
index 837a7a3..2c9ea5f 100644 (file)
@@ -32,7 +32,7 @@ C Virtual-bond lenghts
      & iphi_end,iphid_start,iphid_end,ibond_start,ibond_end,
      & ibondp_start,ibondp_end,ivec_start,ivec_end,iset_start,iset_end,
      & iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start,
-     & iint_end,iphi1_start,iphi1_end,
+     & iint_end,iphi1_start,iphi1_end,itau_start,itau_end,
      & ibond_displ(0:max_fg_procs-1),ibond_count(0:max_fg_procs-1),
      & ithet_displ(0:max_fg_procs-1),ithet_count(0:max_fg_procs-1),
      & iphi_displ(0:max_fg_procs-1),iphi_count(0:max_fg_procs-1),
@@ -46,7 +46,7 @@ C Virtual-bond lenghts
      & ibondp_start,ibondp_end,ivec_start,ivec_end,iset_start,iset_end,
      & iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start,
      & iint_end,iphi1_start,iphi1_end,iint_count,iint_displ,ivec_displ,
-     & ivec_count,iset_displ,
+     & ivec_count,iset_displ,itau_start,itau_end,
      & iset_count,ibond_displ,ibond_count,ithet_displ,ithet_count,
      & iphi_displ,iphi_count,iphi1_displ,iphi1_count
 C Inverses of the actual virtual bond lengths
index a28f621..72a9a1d 100644 (file)
@@ -1,6 +1,19 @@
-C Parameters of the SCCOR term
-      double precision v1sccor,v2sccor
-      integer nterm_sccor
-      common/sccor/v1sccor(maxterm_sccor,20,20),
-     &    v2sccor(maxterm_sccor,20,20),
-     &    nterm_sccor
+cc Parameters of the SCCOR term
+      double precision v1sccor,v2sccor,vlor1sccor,
+     &                 vlor2sccor,vlor3sccor,gloc_sc,
+     &                 dcostau,dsintau,dtauangle,dcosomicron,
+     &                 domicron
+      integer nterm_sccor,isccortyp,nsccortyp,nlor_sccor
+      common/sccor/v1sccor(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp),
+     &    v2sccor(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp),
+     &    v0sccor(maxterm_sccor,-ntyp:ntyp),
+     &    nterm_sccor(-ntyp:ntyp,-ntyp:ntyp),isccortyp(-ntyp:ntyp),
+     &    nsccortyp,
+     &    nlor_sccor(-ntyp:ntyp,-ntyp:ntyp),
+     &    vlor1sccor(maxterm_sccor,20,20),
+     &    vlor2sccor(maxterm_sccor,20,20),
+     &    vlor3sccor(maxterm_sccor,20,20),gloc_sc(3,0:maxres2,10),
+     &    dcostau(3,3,3,maxres2),dsintau(3,3,3,maxres2),
+     &    dtauangle(3,3,3,maxres2),dcosomicron(3,3,3,maxres2),
+     &    domicron(3,3,3,maxres2)
+
index 71158b8..1ab0a16 100644 (file)
@@ -3,10 +3,12 @@ C Store the geometric variables in the following COMMON block.
      &        mask_theta,mask_phi,mask_side
       double precision theta,phi,alph,omeg,varsave,esave,varall,vbld,
      &          thetaref,phiref,costtab,sinttab,cost2tab,sint2tab,
+     &          tauangle,omicron,
      &          xxtab,yytab,zztab,xxref,yyref,zzref
       common /var/ theta(maxres),phi(maxres),alph(maxres),omeg(maxres),
      &          vbld(2*maxres),thetaref(maxres),phiref(maxres),
      &          costtab(maxres), sinttab(maxres), cost2tab(maxres),
+     &          omicron(2,maxres),tauangle(3,maxres),
      &          sint2tab(maxres),xxtab(maxres),yytab(maxres),
      &          zztab(maxres),xxref(maxres),yyref(maxres),zzref(maxres),
      &          ialph(maxres,2),ivar(4*maxres2),ntheta,nphi,nside,nvar
index 3225a09..94d9e45 100644 (file)
@@ -6,7 +6,7 @@
 ********************************************************************************
 C Max. number of processors.
       integer maxprocs
-      parameter (maxprocs=2048)
+      parameter (maxprocs=1024)
 C Max. number of fine-grain processors
       integer max_fg_procs
 c      parameter (max_fg_procs=maxprocs)
@@ -16,7 +16,7 @@ C Max. number of coarse-grain processors
       parameter (max_cg_procs=maxprocs)
 C Max. number of AA residues
       integer maxres
-      parameter (maxres=800)
+      parameter (maxres=600)
 C Appr. max. number of interaction sites
       integer maxres2,maxres6,mmaxres2
       parameter (maxres2=2*maxres,maxres6=6*maxres)
index 3fc2e62..39882d3 100644 (file)
@@ -520,7 +520,20 @@ c-------------------------------------------------------------------------
      &     +(c(j,i+1)-c(j,i))/dnorm2)
         enddo
         be=0.0D0
-        if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
+        if (i.gt.2) then
+        if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
+        if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
+         tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
+        endif
+        if (itype(i-1).ne.10) then
+         tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
+         omicron(1,i)=alpha(i-2,i-1,i-1+nres)
+         omicron(2,i)=alpha(i-1+nres,i-1,i)
+        endif
+        if (itype(i).ne.10) then
+         tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
+        endif
+        endif
         omeg(i)=beta(nres+i,i,maxres2,i+1)
         alph(i)=alpha(nres+i,i,maxres2)
         theta(i+1)=alpha(i-1,i,i+1)
index 015a2b3..73bb4d6 100644 (file)
@@ -5866,28 +5866,52 @@ c        amino-acid residues.
 C Set lprn=.true. for debugging
       lprn=.false.
 c      lprn=.true.
-c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
+c      write (iout,*) "EBACK_SC_COR",itau_start,itau_end
       esccor=0.0D0
-      do i=iphi_start,iphi_end
+      do i=itau_start,itau_end
         esccor_ii=0.0D0
-        itori=itype(i-2)
-        itori1=itype(i-1)
+        isccori=isccortyp(itype(i-2))
+        isccori1=isccortyp(itype(i-1))
+c      write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1)
         phii=phi(i)
+        do intertyp=1,3 !intertyp
+cc Added 09 May 2012 (Adasko)
+cc  Intertyp means interaction type of backbone mainchain correlation: 
+c   1 = SC...Ca...Ca...Ca
+c   2 = Ca...Ca...Ca...SC
+c   3 = SC...Ca...Ca...SCi
         gloci=0.0D0
-        do j=1,nterm_sccor
-          v1ij=v1sccor(j,itori,itori1)
-          v2ij=v2sccor(j,itori,itori1)
-          cosphi=dcos(j*phii)
-          sinphi=dsin(j*phii)
+        if (((intertyp.eq.3).and.((itype(i-2).eq.10).or.
+     &      (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or.
+     &      (itype(i-1).eq.ntyp1)))
+     &    .or. ((intertyp.eq.1).and.((itype(i-2).eq.10)
+     &     .or.(itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)
+     &     .or.(itype(i).eq.ntyp1)))
+     &    .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or.
+     &      (itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+     &      (itype(i-3).eq.ntyp1)))) cycle
+        if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle
+        if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1))
+     & cycle
+       do j=1,nterm_sccor(isccori,isccori1)
+          v1ij=v1sccor(j,intertyp,isccori,isccori1)
+          v2ij=v2sccor(j,intertyp,isccori,isccori1)
+          cosphi=dcos(j*tauangle(intertyp,i))
+          sinphi=dsin(j*tauangle(intertyp,i))
           esccor=esccor+v1ij*cosphi+v2ij*sinphi
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
+c      write (iout,*) "EBACK_SC_COR",i,esccor,intertyp
+        gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
         if (lprn)
      &  write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
-     &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
-     &  (v1sccor(j,itori,itori1),j=1,6),(v2sccor(j,itori,itori1),j=1,6)
+     &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,isccori,isccori1,
+     &  (v1sccor(j,intertyp,isccori,isccori1),j=1,6)
+     & ,(v2sccor(j,intertyp,isccori,isccori1),j=1,6)
         gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
+       enddo !intertyp
       enddo
+
       return
       end
 c----------------------------------------------------------------------------
index 25d1b12..569120f 100644 (file)
@@ -337,6 +337,7 @@ C-------------------------------------------------------------------------
       include 'COMMON.CHAIN'
       include 'COMMON.VAR'
       include 'COMMON.MD_'
+      include 'COMMON.SCCOR'
 C
 C Initialize Cartesian-coordinate gradient
 C
@@ -374,6 +375,9 @@ C
           gradx(j,i,icg)=0.0d0
           gscloc(j,i)=0.0d0
           gsclocx(j,i)=0.0d0
+          do intertyp=1,3
+           gloc_sc(intertyp,i,icg)=0.0d0
+          enddo
         enddo
       enddo
 C
index b3affab..d2d8f9a 100644 (file)
@@ -565,6 +565,9 @@ C Partition local interactions
       iphi_end=iturn3_end+2
       iturn3_start=iturn3_start-1
       iturn3_end=iturn3_end-1
+      call int_bounds(nres-3,itau_start,itau_end)
+      itau_start=itau_start+3
+      itau_end=itau_end+3
       call int_bounds(nres-3,iphi1_start,iphi1_end)
       iphi1_start=iphi1_start+3
       iphi1_end=iphi1_end+3
@@ -1091,6 +1094,8 @@ c        write (iout,*) "MPI_ROTAT2",MPI_ROTAT2
       idihconstr_end=ndih_constr
       iphid_start=iphi_start
       iphid_end=iphi_end-1
+      itau_start=4
+      itau_end=nres
       ibond_start=2
       ibond_end=nres-1
       ibondp_start=nnt+1
index 97324ec..18f7d54 100644 (file)
@@ -13,9 +13,9 @@ c-------------------------------------------------------------
       include 'COMMON.INTERACT'
       include 'COMMON.MD_'
       include 'COMMON.IOUNITS'
-      
+      include 'COMMON.SCCOR'       
 c   calculating dE/ddc1      
-       if (nres.lt.3) return
+       if (nres.lt.3) go to 18
        do j=1,3
          gcart(j,1)=gcart(j,1)+gloc(1,icg)*dphi(j,1,4)
      &     +gloc(nres-2,icg)*dtheta(j,1,3)      
@@ -113,6 +113,160 @@ c   The side-chain vector derivatives
             enddo
          endif     
        enddo                                                                                                                                                   
+c----------------------------------------------------------------------
+C INTERTYP=1 SC...Ca...Ca...Ca
+C INTERTYP=2 Ca...Ca...Ca...SC
+C INTERTYP=3 SC...Ca...Ca...SC
+c   calculating dE/ddc1      
+  18   continue
+c       do i=1,nres
+c       gloc(i,icg)=0.0D0
+c          write (iout,*) "poczotkoawy",i,gloc_sc(1,i,icg)
+c       enddo
+       if (nres.lt.2) return
+       if ((nres.lt.3).and.(itype(1).eq.10)) return
+       if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then
+        do j=1,3
+cc Derviative was calculated for oposite vector of side chain therefore
+c there is "-" sign before gloc_sc
+         gxcart(j,1)=gxcart(j,1)-gloc_sc(1,0,icg)*
+     &     dtauangle(j,1,1,3)
+         gcart(j,1)=gcart(j,1)+gloc_sc(1,0,icg)*
+     &     dtauangle(j,1,2,3)
+          if ((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then
+         gxcart(j,1)= gxcart(j,1)
+     &               -gloc_sc(3,0,icg)*dtauangle(j,3,1,3)
+         gcart(j,1)=gcart(j,1)+gloc_sc(3,0,icg)*
+     &       dtauangle(j,3,2,3)
+          endif
+       enddo
+       endif
+         if ((nres.ge.3).and.(itype(3).ne.10).and.(itype(3).ne.ntyp1))
+     & then
+         do j=1,3
+         gcart(j,1)=gcart(j,1)+gloc_sc(2,1,icg)*dtauangle(j,2,1,4)
+         enddo
+         endif
+c   As potetnial DO NOT depend on omicron anlge their derivative is
+c   ommited 
+c     &     +gloc_sc(intertyp,nres-2,icg)*dtheta(j,1,3)  
+
+c     Calculating the remainder of dE/ddc2
+       do j=1,3
+         if((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then
+           if (itype(1).ne.10) gxcart(j,2)=gxcart(j,2)+
+     &                         gloc_sc(3,0,icg)*dtauangle(j,3,3,3)
+        if ((itype(3).ne.10).and.(nres.ge.3).and.(itype(3).ne.ntyp1))
+     &   then
+           gxcart(j,2)=gxcart(j,2)-gloc_sc(3,1,icg)*dtauangle(j,3,1,4)
+cc                  the   - above is due to different vector direction
+           gcart(j,2)=gcart(j,2)+gloc_sc(3,1,icg)*dtauangle(j,3,2,4)
+          endif
+          if (nres.gt.3) then
+           gxcart(j,2)=gxcart(j,2)-gloc_sc(1,1,icg)*dtauangle(j,1,1,4)
+cc                  the   - above is due to different vector direction
+           gcart(j,2)=gcart(j,2)+gloc_sc(1,1,icg)*dtauangle(j,1,2,4)
+c          write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,2,4),"gcart"
+c           write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx"
+          endif
+         endif
+         if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then
+          gcart(j,2)=gcart(j,2)+gloc_sc(1,0,icg)*dtauangle(j,1,3,3)
+c           write(iout,*)  gloc_sc(1,0,icg),dtauangle(j,1,3,3)
+        endif
+         if ((itype(3).ne.10).and.(nres.ge.3)) then
+          gcart(j,2)=gcart(j,2)+gloc_sc(2,1,icg)*dtauangle(j,2,2,4)
+c           write(iout,*) gloc_sc(2,1,icg),dtauangle(j,2,2,4)
+         endif
+         if ((itype(4).ne.10).and.(nres.ge.4)) then
+          gcart(j,2)=gcart(j,2)+gloc_sc(2,2,icg)*dtauangle(j,2,1,5)
+c           write(iout,*) gloc_sc(2,2,icg),dtauangle(j,2,1,5)
+         endif
+
+c      write(iout,*) gcart(j,2),itype(2),itype(1),itype(3), "gcart2"
+       enddo
+c    If there are more than five residues
+      if(nres.ge.5) then                        
+        do i=3,nres-2
+         do j=1,3
+c          write(iout,*) "before", gcart(j,i)
+          if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+          gxcart(j,i)=gxcart(j,i)+gloc_sc(2,i-2,icg)
+     &    *dtauangle(j,2,3,i+1)
+     &    -gloc_sc(1,i-1,icg)*dtauangle(j,1,1,i+2)
+          gcart(j,i)=gcart(j,i)+gloc_sc(1,i-1,icg)
+     &    *dtauangle(j,1,2,i+2)
+c                   write(iout,*) "new",j,i,
+c     &  gcart(j,i),gloc_sc(1,i-1,icg),dtauangle(j,1,2,i+2)
+          if (itype(i-1).ne.10) then
+           gxcart(j,i)=gxcart(j,i)+gloc_sc(3,i-2,icg)
+     &*dtauangle(j,3,3,i+1)
+          endif
+          if (itype(i+1).ne.10) then
+           gxcart(j,i)=gxcart(j,i)-gloc_sc(3,i-1,icg)
+     &*dtauangle(j,3,1,i+2)
+           gcart(j,i)=gcart(j,i)+gloc_sc(3,i-1,icg)
+     &*dtauangle(j,3,2,i+2)
+          endif
+          endif
+          if (itype(i-1).ne.10) then
+           gcart(j,i)=gcart(j,i)+gloc_sc(1,i-2,icg)*
+     &     dtauangle(j,1,3,i+1)
+          endif
+          if (itype(i+1).ne.10) then
+           gcart(j,i)=gcart(j,i)+gloc_sc(2,i-1,icg)*
+     &     dtauangle(j,2,2,i+2)
+c          write(iout,*) "numer",i,gloc_sc(2,i-1,icg),
+c     &    dtauangle(j,2,2,i+2)
+          endif
+          if (itype(i+2).ne.10) then
+           gcart(j,i)=gcart(j,i)+gloc_sc(2,i,icg)*
+     &     dtauangle(j,2,1,i+3)
+          endif
+         enddo
+        enddo
+      endif     
+c  Setting dE/ddnres-1       
+      if(nres.ge.4) then
+         do j=1,3
+         if ((itype(nres-1).ne.10).and.(itype(nres-1).ne.ntyp1)) then
+         gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(2,nres-3,icg)
+     &    *dtauangle(j,2,3,nres)
+c          write (iout,*) "gxcart(nres-1)", gloc_sc(2,nres-3,icg),
+c     &     dtauangle(j,2,3,nres), gxcart(j,nres-1)
+         if (itype(nres-2).ne.10) then
+        gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg)
+     &    *dtauangle(j,3,3,nres)
+          endif
+         if ((itype(nres).ne.10).and.(itype(nres).ne.ntyp1)) then
+        gxcart(j,nres-1)=gxcart(j,nres-1)-gloc_sc(3,nres-2,icg)
+     &    *dtauangle(j,3,1,nres+1)
+        gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(3,nres-2,icg)
+     &    *dtauangle(j,3,2,nres+1)
+          endif
+         endif
+         if ((itype(nres-2).ne.10).and.(itype(nres-2).ne.ntyp1)) then
+            gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(1,nres-3,icg)*
+     &   dtauangle(j,1,3,nres)
+         endif
+          if ((itype(nres).ne.10).and.(itype(nres).ne.ntyp1)) then
+            gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(2,nres-2,icg)*
+     &     dtauangle(j,2,2,nres+1)
+c           write (iout,*) "gcart(nres-1)", gloc_sc(2,nres-2,icg),
+c     &     dtauangle(j,2,2,nres+1), itype(nres-1),itype(nres)
+           endif
+         enddo
+      endif
+c  Settind dE/ddnres       
+       if ((nres.ge.3).and.(itype(nres).ne.10).and.
+     &    (itype(nres).ne.ntyp1))then
+       do j=1,3
+        gxcart(j,nres)=gxcart(j,nres)+gloc_sc(3,nres-2,icg)
+     & *dtauangle(j,3,3,nres+1)+gloc_sc(2,nres-2,icg)
+     & *dtauangle(j,2,3,nres+1)
+        enddo
+       endif
+c   The side-chain vector derivatives
       return
       end      
        
index 5fea875..9eb0098 100644 (file)
@@ -12,6 +12,7 @@
       include 'COMMON.DERIV'
       include 'COMMON.IOUNITS'
       include 'COMMON.LOCAL'
+      include 'COMMON.SCCOR'
       double precision dcostheta(3,2,maxres),
      & dcosphi(3,3,maxres),dsinphi(3,3,maxres),
      & dcosalpha(3,3,maxres),dcosomega(3,3,maxres),
@@ -44,7 +45,41 @@ c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
           dtheta(j,2,i)=-1/sint*dcostheta(j,2,i)     
         enddo
       enddo
-      
+#if defined(MPI) && defined(PARINTDER)
+c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
+      do i=max0(ithet_start-1,3),ithet_end
+#else
+      do i=3,nres
+#endif
+      if ((itype(i-1).ne.10).and.(itype(i-1).ne.ntyp1)) then
+        cost1=dcos(omicron(1,i))
+        sint1=sqrt(1-cost1*cost1)
+        cost2=dcos(omicron(2,i))
+        sint2=sqrt(1-cost2*cost2)
+       do j=1,3
+CC Calculate derivative over first omicron (Cai-2,Cai-1,SCi-1) 
+          dcosomicron(j,1,1,i)=-(dc_norm(j,i-1+nres)+
+     &    cost1*dc_norm(j,i-2))/
+     &    vbld(i-1)
+          domicron(j,1,1,i)=-1/sint1*dcosomicron(j,1,1,i)
+          dcosomicron(j,1,2,i)=-(dc_norm(j,i-2)
+     &    +cost1*(dc_norm(j,i-1+nres)))/
+     &    vbld(i-1+nres)
+          domicron(j,1,2,i)=-1/sint1*dcosomicron(j,1,2,i)
+CC Calculate derivative over second omicron Sci-1,Cai-1 Cai
+CC Looks messy but better than if in loop
+          dcosomicron(j,2,1,i)=-(-dc_norm(j,i-1+nres)
+     &    +cost2*dc_norm(j,i-1))/
+     &    vbld(i)
+          domicron(j,2,1,i)=-1/sint2*dcosomicron(j,2,1,i)
+          dcosomicron(j,2,2,i)=-(dc_norm(j,i-1)
+     &     +cost2*(-dc_norm(j,i-1+nres)))/
+     &    vbld(i-1+nres)
+c          write(iout,*) "vbld", i,itype(i),vbld(i-1+nres)
+          domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i)
+        enddo
+       endif
+      enddo      
 c Derivatives of phi:
 c If phi is 0 or 180 degrees, then the formulas 
 c have to be derived by power series expansion of the
@@ -109,6 +144,227 @@ c   Obtaining the gamma derivatives from cosine derivative
          enddo
         endif                                                                                           
       enddo
+Calculate derivative of Tauangle
+#ifdef PARINTDER
+      do i=itau_start,itau_end
+#else
+      do i=3,nres
+#endif
+       if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle
+c       if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10).or.
+c     &     (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) cycle
+cc dtauangle(j,intertyp,dervityp,residue number)
+cc INTERTYP=1 SC...Ca...Ca..Ca
+c the conventional case
+        sint=dsin(theta(i))
+        sint1=dsin(omicron(2,i-1))
+        sing=dsin(tauangle(1,i))
+        cost=dcos(theta(i))
+        cost1=dcos(omicron(2,i-1))
+        cosg=dcos(tauangle(1,i))
+        do j=1,3
+        dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
+cc       write(iout,*) dc_norm2(j,i-2+nres),"dcnorm"
+        enddo
+        scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1))
+        fac0=1.0d0/(sint1*sint)
+        fac1=cost*fac0
+        fac2=cost1*fac0
+        fac3=cosg*cost1/(sint1*sint1)
+        fac4=cosg*cost/(sint*sint)
+cc         write(iout,*) "faki",fac0,fac1,fac2,fac3,fac4
+c    Obtaining the gamma derivatives from sine derivative                                
+       if (tauangle(1,i).gt.-pi4.and.tauangle(1,i).le.pi4.or.
+     &     tauangle(1,i).gt.pi34.and.tauangle(1,i).le.pi.or.
+     &     tauangle(1,i).gt.-pi.and.tauangle(1,i).le.-pi34) then
+         call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
+         call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1),vp2)
+         call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
+        do j=1,3
+            ctgt=cost/sint
+            ctgt1=cost1/sint1
+            cosg_inv=1.0d0/cosg
+            dsintau(j,1,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
+     &-(fac0*vp1(j)+sing*(dc_norm2(j,i-2+nres)))
+     & *vbld_inv(i-2+nres)
+            dtauangle(j,1,1,i)=cosg_inv*dsintau(j,1,1,i)
+            dsintau(j,1,2,i)=
+     &        -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*dtheta(j,1,i))
+     &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+c            write(iout,*) "dsintau", dsintau(j,1,2,i)
+            dtauangle(j,1,2,i)=cosg_inv*dsintau(j,1,2,i)
+c Bug fixed 3/24/05 (AL)
+            dsintau(j,1,3,i)=-sing*ctgt*dtheta(j,2,i)
+     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
+c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
+            dtauangle(j,1,3,i)=cosg_inv*dsintau(j,1,3,i)
+         enddo
+c   Obtaining the gamma derivatives from cosine derivative
+        else
+           do j=1,3
+           dcostau(j,1,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
+     &     dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1)-scalp*
+     &     (dc_norm2(j,i-2+nres)))/vbld(i-2+nres)
+           dtauangle(j,1,1,i)=-1/sing*dcostau(j,1,1,i)
+           dcostau(j,1,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
+     &     dcostheta(j,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
+     &     dcostheta(j,1,i)
+           dtauangle(j,1,2,i)=-1/sing*dcostau(j,1,2,i)
+           dcostau(j,1,3,i)=fac2*dcostheta(j,2,i)+fac4*
+     &     dcostheta(j,2,i)-fac0*(-dc_norm(j,i-2+nres)-scalp*
+     &     dc_norm(j,i-1))/vbld(i)
+           dtauangle(j,1,3,i)=-1/sing*dcostau(j,1,3,i)
+c         write (iout,*) "else",i
+         enddo
+        endif
+c        do k=1,3                 
+c        write(iout,*) "tu",i,k,(dtauangle(j,1,k,i),j=1,3)        
+c        enddo                
+      enddo
+CC Second case Ca...Ca...Ca...SC
+#ifdef PARINTDER
+      do i=itau_start,itau_end
+#else
+      do i=4,nres
+#endif
+       if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or.
+     &    (itype(i-2).eq.ntyp1).or.(itype(i-3).eq.ntyp1)) cycle
+c the conventional case
+        sint=dsin(omicron(1,i))
+        sint1=dsin(theta(i-1))
+        sing=dsin(tauangle(2,i))
+        cost=dcos(omicron(1,i))
+        cost1=dcos(theta(i-1))
+        cosg=dcos(tauangle(2,i))
+c        do j=1,3
+c        dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
+c        enddo
+        scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1+nres))
+        fac0=1.0d0/(sint1*sint)
+        fac1=cost*fac0
+        fac2=cost1*fac0
+        fac3=cosg*cost1/(sint1*sint1)
+        fac4=cosg*cost/(sint*sint)
+c    Obtaining the gamma derivatives from sine derivative                                
+       if (tauangle(2,i).gt.-pi4.and.tauangle(2,i).le.pi4.or.
+     &     tauangle(2,i).gt.pi34.and.tauangle(2,i).le.pi.or.
+     &     tauangle(2,i).gt.-pi.and.tauangle(2,i).le.-pi34) then
+         call vecpr(dc_norm2(1,i-1+nres),dc_norm(1,i-2),vp1)
+         call vecpr(dc_norm(1,i-3),dc_norm(1,i-1+nres),vp2)
+         call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
+        do j=1,3
+            ctgt=cost/sint
+            ctgt1=cost1/sint1
+            cosg_inv=1.0d0/cosg
+            dsintau(j,2,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
+     &        +(fac0*vp1(j)-sing*dc_norm(j,i-3))*vbld_inv(i-2)
+c       write(iout,*) i,j,dsintau(j,2,1,i),sing*ctgt1*dtheta(j,1,i-1),
+c     &fac0*vp1(j),sing*dc_norm(j,i-3),vbld_inv(i-2),"dsintau(2,1)"
+            dtauangle(j,2,1,i)=cosg_inv*dsintau(j,2,1,i)
+            dsintau(j,2,2,i)=
+     &        -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*domicron(j,1,1,i))
+     &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+c            write(iout,*) "sprawdzenie",i,j,sing*ctgt1*dtheta(j,2,i-1),
+c     & sing*ctgt*domicron(j,1,2,i),
+c     & (fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+            dtauangle(j,2,2,i)=cosg_inv*dsintau(j,2,2,i)
+c Bug fixed 3/24/05 (AL)
+            dsintau(j,2,3,i)=-sing*ctgt*domicron(j,1,2,i)
+     &       +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))*vbld_inv(i-1+nres)
+c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
+            dtauangle(j,2,3,i)=cosg_inv*dsintau(j,2,3,i)
+         enddo
+c   Obtaining the gamma derivatives from cosine derivative
+        else
+           do j=1,3
+           dcostau(j,2,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
+     &     dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
+     &     dc_norm(j,i-3))/vbld(i-2)
+           dtauangle(j,2,1,i)=-1/sing*dcostau(j,2,1,i)
+           dcostau(j,2,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
+     &     dcosomicron(j,1,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
+     &     dcosomicron(j,1,1,i)
+           dtauangle(j,2,2,i)=-1/sing*dcostau(j,2,2,i)
+           dcostau(j,2,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
+     &     dcosomicron(j,1,2,i)-fac0*(dc_norm(j,i-3)-scalp*
+     &     dc_norm(j,i-1+nres))/vbld(i-1+nres)
+           dtauangle(j,2,3,i)=-1/sing*dcostau(j,2,3,i)
+c        write(iout,*) i,j,"else", dtauangle(j,2,3,i) 
+         enddo
+        endif                                    
+      enddo
+
+CCC third case SC...Ca...Ca...SC
+#ifdef PARINTDER
+
+      do i=itau_start,itau_end
+#else
+      do i=3,nres
+#endif
+c the conventional case
+      if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or.
+     &(itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle
+        sint=dsin(omicron(1,i))
+        sint1=dsin(omicron(2,i-1))
+        sing=dsin(tauangle(3,i))
+        cost=dcos(omicron(1,i))
+        cost1=dcos(omicron(2,i-1))
+        cosg=dcos(tauangle(3,i))
+        do j=1,3
+        dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
+c        dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
+        enddo
+        scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres))
+        fac0=1.0d0/(sint1*sint)
+        fac1=cost*fac0
+        fac2=cost1*fac0
+        fac3=cosg*cost1/(sint1*sint1)
+        fac4=cosg*cost/(sint*sint)
+c    Obtaining the gamma derivatives from sine derivative                                
+       if (tauangle(3,i).gt.-pi4.and.tauangle(3,i).le.pi4.or.
+     &     tauangle(3,i).gt.pi34.and.tauangle(3,i).le.pi.or.
+     &     tauangle(3,i).gt.-pi.and.tauangle(3,i).le.-pi34) then
+         call vecpr(dc_norm(1,i-1+nres),dc_norm(1,i-2),vp1)
+         call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres),vp2)
+         call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
+        do j=1,3
+            ctgt=cost/sint
+            ctgt1=cost1/sint1
+            cosg_inv=1.0d0/cosg
+            dsintau(j,3,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
+     &        -(fac0*vp1(j)-sing*dc_norm(j,i-2+nres))
+     &        *vbld_inv(i-2+nres)
+            dtauangle(j,3,1,i)=cosg_inv*dsintau(j,3,1,i)
+            dsintau(j,3,2,i)=
+     &        -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*domicron(j,1,1,i))
+     &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+            dtauangle(j,3,2,i)=cosg_inv*dsintau(j,3,2,i)
+c Bug fixed 3/24/05 (AL)
+            dsintau(j,3,3,i)=-sing*ctgt*domicron(j,1,2,i)
+     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))
+     &        *vbld_inv(i-1+nres)
+c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
+            dtauangle(j,3,3,i)=cosg_inv*dsintau(j,3,3,i)
+         enddo
+c   Obtaining the gamma derivatives from cosine derivative
+        else
+           do j=1,3
+           dcostau(j,3,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
+     &     dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
+     &     dc_norm2(j,i-2+nres))/vbld(i-2+nres)
+           dtauangle(j,3,1,i)=-1/sing*dcostau(j,3,1,i)
+           dcostau(j,3,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
+     &     dcosomicron(j,1,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
+     &     dcosomicron(j,1,1,i)
+           dtauangle(j,3,2,i)=-1/sing*dcostau(j,3,2,i)
+           dcostau(j,3,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
+     &     dcosomicron(j,1,2,i)-fac0*(dc_norm2(j,i-2+nres)-scalp*
+     &     dc_norm(j,i-1+nres))/vbld(i-1+nres)
+           dtauangle(j,3,3,i)=-1/sing*dcostau(j,3,3,i)
+c          write(iout,*) "else",i 
+         enddo
+        endif                                                                                            
+      enddo
 #ifdef CRYST_SC
 c   Derivatives of side-chain angles alpha and omega
 #if defined(MPI) && defined(PARINTDER)
index 0a99250..3955c0c 100644 (file)
@@ -553,28 +553,133 @@ C
       enddo
       endif
 #endif
+C Read of Side-chain backbone correlation parameters
+C Modified 11 May 2012 by Adasko
+CCC
 C
 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
 C         correlation energies.
 C
-      read (isccor,*,end=119,err=119) nterm_sccor
-      do i=1,20
-       do j=1,20
-          read (isccor,'(a)')
-         do k=1,nterm_sccor
-           read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j),
-     &        v2sccor(k,i,j) 
+      read (isccor,*,end=119,err=119) nsccortyp
+#ifdef SCCORPDB
+      read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
+C For D-aminoacid uncomment
+C      do i=-ntyp,-1
+C        isccortyp(i)=-isccortyp(-i)
+C      enddo
+      iscprol=isccortyp(20)
+c      write (iout,*) 'ntortyp',ntortyp
+      maxinter=3
+cc maxinter is maximum interaction sites
+      do l=1,maxinter
+      do i=1,nsccortyp
+        do j=1,nsccortyp
+          read (isccor,*,end=119,err=119)
+     &nterm_sccor(i,j),nlor_sccor(i,j)
+          v0ijsccor=0.0d0
+          si=-1.0d0
+C          nterm_sccor(-i,j)=nterm_sccor(i,j)
+C          nterm_sccor(-i,-j)=nterm_sccor(i,j)
+C          nterm_sccor(i,-j)=nterm_sccor(i,j)
+          do k=1,nterm_sccor(i,j)
+            read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
+     &    ,v2sccor(k,l,i,j)
+C            if (j.eq.iscprol) then
+C             if (i.eq.isccortyp(10)) then
+C             v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
+C             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+C             else
+C             v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
+C     &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
+C             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
+C     &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
+C             v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
+C             v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
+C             v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
+C             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
+C             endif
+CC           else
+C             if (i.eq.isccortyp(10)) then
+C             v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
+C             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+C             else
+C               if (j.eq.isccortyp(10)) then
+C             v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
+C             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
+C               else
+C             v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
+C             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+C             v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
+C             v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
+C             v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
+C             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
+C                endif
+C               endif
+C            endif
+            v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
+            si=-si
+          enddo
+          do k=1,nlor_sccor(i,j)
+            read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
+     &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
+            v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
+     &(1+vlor3sccor(k,i,j)**2)
+          enddo
+          v0sccor(i,j)=v0ijsccor
+        enddo
+      enddo
+      enddo
+      close (isccor)
+#else
+      read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
+c      write (iout,*) 'ntortyp',ntortyp
+      maxinter=3
+cc maxinter is maximum interaction sites
+      do l=1,maxinter
+      do i=1,nsccortyp
+        do j=1,nsccortyp
+          read (isccor,*,end=113,err=113)
+     & nterm_sccor(i,j),nlor_sccor(i,j)
+          v0ijsccor=0.0d0
+          si=-1.0d0
+
+          do k=1,nterm_sccor(i,j)
+            read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
+     &    ,v2sccor(k,l,i,j)
+            v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
+            si=-si
           enddo
+          do k=1,nlor_sccor(i,j)
+            read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
+     &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
+            v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
+     &(1+vlor3sccor(k,i,j)**2)
+          enddo
+          v0sccor(i,j)=v0ijsccor
         enddo
       enddo
+      enddo
       close (isccor)
+
+#endif      
       if (lprint) then
-       write (iout,'(/a/)') 'Torsional constants of SCCORR:'
-       do i=1,20
-         do j=1,20
+       write (iout,'(/a/)') 'Torsional constants of SCCORR:'
+C       do i=1,20
+C         do j=1,20
+C        write (iout,'(/a/)') 'Torsional constants:'
+        do i=1,nsccortyp
+          do j=1,nsccortyp
             write (iout,*) 'ityp',i,' jtyp',j
-            do k=1,nterm_sccor
-             write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
+C            do k=1,nterm_sccor
+C             write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
+C            write (iout,*) 'Fourier constants'
+            do k=1,nterm_sccor(i,j)
+      write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
+            enddo
+            write (iout,*) 'Lorenz constants'
+            do k=1,nlor_sccor(i,j)
+              write (iout,'(3(1pe15.5))')
+     &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
             enddo
           enddo
         enddo
index eb4ba3f..8dec8a8 100644 (file)
@@ -251,6 +251,7 @@ ctest          stop
         vbld_inv(i+1)=1.0d0/vbld(i+1)
         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
+
       enddo
 c      if (unres_pdb) then
 c        if (itype(1).eq.21) then
index e178afe..12ebe53 100644 (file)
@@ -78,7 +78,7 @@ c    If there are more than five residues
         enddo
       endif    
 c  Setting dE/ddnres-2      
-      write(iout,*) "ATUCHUJ?" 
+C      write(iout,*) "ATUCHUJ?" 
       if(nres.gt.5) then
          do j=1,3
            gcart(j,nres-2)=gcart(j,nres-2)+gloc(nres-4,icg)*
@@ -116,7 +116,7 @@ c   The side-chain vector derivatives
             enddo
          endif     
        enddo
-       write(iout,*) "TU DOCHODZE"
+C       write(iout,*) "TU DOCHODZE"
                                                                                                                                                        
 c----------------------------------------------------------------------
 C INTERTYP=1 SC...Ca...Ca...Ca
index 4c9a018..d1ddcd6 100644 (file)
      & v1sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm),
      & v2sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
       integer nlob_all(ntyp1,max_parm),
-     & nlor_all(-maxtor:maxtor,-maxtor:maxtor,2,max_parm),
-     & nterm_all(-maxtor:maxtor,-maxtor:maxtor,2,max_parm),
+     & nlor_all(-maxtor:maxtor,-maxtor:maxtor,max_parm),
+     & nterm_all(-maxtor:maxtor,-maxtor:maxtor,max_parm),
      & ntermd1_all(-maxtor:maxtor,-maxtor:maxtor,
-     & -maxtor:maxtor,2,max_parm),
+     & -maxtor:maxtor,max_parm),
      & ntermd2_all(-maxtor:maxtor,-maxtor:maxtor,
-     & -maxtor:maxtor,2,max_parm),
+     & -maxtor:maxtor,max_parm),
      & nbondterm_all(ntyp,max_parm),nthetyp_all(max_parm),
      & ithetyp_all(ntyp1,max_parm),ntheterm_all(max_parm),
      & ntheterm2_all(max_parm),ntheterm3_all(max_parm),
index 42ace37..a248d99 100644 (file)
@@ -2,7 +2,7 @@
      & sigc0,dsc,dsc_inv,bsc,censc,gaussc,dsc0,vbl,vblinv,vblinv2,
      & vbl_cis,vbl0,vbld_inv
       integer nlob,loc_start,loc_end,ithet_start,ithet_end,
-     & iphi_start,iphi_end
+     & iphi_start,iphi_end,itau_start,itau_end
 C Parameters of the virtual-bond-angle probability distribution
       common /thetas/ a0thet(ntyp),athet(2,ntyp),bthet(2,ntyp),
      &  polthet(0:3,ntyp),gthet(3,ntyp),theta0(ntyp),sig0(ntyp),
@@ -29,17 +29,8 @@ C Parameters of the side-chain probability distribution
      &  censc(3,maxlob,ntyp),gaussc(3,3,maxlob,ntyp),dsc0(ntyp1),
      &    nlob(ntyp1)
 C Virtual-bond lenghts
-      double precision vbl,vblinv,vblinv2,vbl_cis,vbl0,vbld_inv
-      integer loc_start,loc_end,ithet_start,ithet_end,iphi_start,
-     & iphi_end,iphid_start,iphid_end,ibond_start,ibond_end,
-     & ibondp_start,ibondp_end,ivec_start,ivec_end,iset_start,iset_end,
-     & iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start,
-     & iint_end,iphi1_start,iphi1_end,itau_start,itau_end
       common /peptbond/ vbl,vblinv,vblinv2,vbl_cis,vbl0
       common /indices/ loc_start,loc_end,ithet_start,ithet_end,
-     & iphi_start,iphi_end,iphid_start,iphid_end,ibond_start,ibond_end,
-     & ibondp_start,ibondp_end,ivec_start,ivec_end,iset_start,iset_end,
-     & iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start,
-     & iint_end,iphi1_start,iphi1_end,itau_start,itau_end
+     &                 iphi_start,iphi_end,itau_start,itau_end
 C Inverses of the actual virtual bond lengths
       common /invlen/ vbld_inv(maxres2)