Dzialajacy INTERTYP=1 i INTERTYP=2
authorAdam Kazimierz Sieradzan <adasko@sun1.chem.univ.gda.pl>
Thu, 14 Jun 2012 11:03:08 +0000 (07:03 -0400)
committerAdam Kazimierz Sieradzan <adasko@sun1.chem.univ.gda.pl>
Thu, 14 Jun 2012 11:03:08 +0000 (07:03 -0400)
14 files changed:
bin/unres/MD/unres_ifort_MPICH_GAB.exe
source/unres/src_MD/COMMON.CHAIN
source/unres/src_MD/COMMON.SCCOR
source/unres/src_MD/COMMON.VAR
source/unres/src_MD/DIMENSIONS
source/unres/src_MD/Makefile_MPICH_ifort
source/unres/src_MD/checkder_p.F
source/unres/src_MD/cinfo.f
source/unres/src_MD/energy_p_new_barrier.F
source/unres/src_MD/gradient_p.F
source/unres/src_MD/int_to_cart.f
source/unres/src_MD/intcartderiv.F
source/unres/src_MD/parmread.F
source/unres/src_MD/tmptmp

index f772e87..e333657 100755 (executable)
Binary files a/bin/unres/MD/unres_ifort_MPICH_GAB.exe and b/bin/unres/MD/unres_ifort_MPICH_GAB.exe differ
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 1c725d7..ca6210f 100644 (file)
@@ -1,11 +1,16 @@
- Parameters of the SCCOR term
+cc Parameters of the SCCOR term
       double precision v1sccor,v2sccor,vlor1sccor,
-     &                 vlor2sccor,vlor3sccor
+     &                 vlor2sccor,vlor3sccor,gloc_sc,
+     &                 dcostau,dsintau,dtauangle,dcosomicron,
+     &                 domicron
       integer nterm_sccor,isccortyp,nsccortyp,nlor_sccor
-      common/torsion/v1sccor(maxterm_sccor,20,20),
-     &    v2sccor(maxterm_sccor,20,20),
+      common/sccor/v1sccor(maxterm_sccor,3,20,20),
+     &    v2sccor(maxterm_sccor,3,20,20),
+     &    v0sccor(maxterm_sccor,20),
      &    nterm_sccor(ntyp,ntyp),isccortyp(ntyp),nsccortyp,
      &    nlor_sccor(ntyp,ntyp),vlor1sccor(maxterm_sccor,20,20),
      &    vlor2sccor(maxterm_sccor,20,20),
-     &    vlor3sccor(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 c2a0a25..d560c87 100644 (file)
@@ -5,7 +5,7 @@ C Store the geometric variables in the following COMMON block.
      &          thetaref,phiref,costtab,sinttab,cost2tab,sint2tab,
      &          xxtab,yytab,zztab,xxref,yyref,zzref
       common /var/ theta(maxres),phi(maxres),alph(maxres),omeg(maxres),
-     &          omicron(2,maxres),tauangle(3,maxres)
+     &          omicron(2,maxres),tauangle(3,maxres),
      &          vbld(2*maxres),thetaref(maxres),phiref(maxres),
      &          costtab(maxres), sinttab(maxres), cost2tab(maxres),
      &          sint2tab(maxres),xxtab(maxres),yytab(maxres),
index 224dade..5151ff7 100644 (file)
@@ -54,7 +54,7 @@ C virtual-bond angle bending potentials
      & mmaxtheterm=maxtheterm)
 c Max number of torsional terms in SCCOR
       integer maxterm_sccor
-      parameter (maxterm_sccor=3)
+      parameter (maxterm_sccor=6)
 C Max. number of lobes in SC distribution
       integer maxlob
       parameter (maxlob=4)
index e33762c..8846882 100644 (file)
@@ -15,7 +15,7 @@ INSTALL_DIR = /users/software/mpich-1.2.7p1_intel-10.1_em64_ssh
 
 FC= ifort
 
-OPT =  -O3 -ip -w 
+OPT =  -g -ip -w -CB 
 
 FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include 
 FFLAGS1 = -c -w -g -d2 -CA -CB -I$(INSTALL_DIR)/include 
index 719770e..4d0379e 100644 (file)
@@ -8,6 +8,7 @@ C Check the gradient of Cartesian coordinates in internal coordinates.
       include 'COMMON.GEO'
       include 'COMMON.LOCAL'
       include 'COMMON.DERIV'
+      include 'COMMON.SCCOR'
       dimension temp(6,maxres),xx(3),gg(3)
       indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
 *
@@ -180,6 +181,7 @@ C Check the gradient of the energy in Cartesian coordinates.
       include 'COMMON.IOUNITS'
       include 'COMMON.VAR'
       include 'COMMON.CONTACTS'
+      include 'COMMON.SCCOR'
       common /srutu/ icall
       dimension ggg(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),g(maxvar)
       dimension grad_s(6,maxres)
@@ -261,6 +263,7 @@ C Check the gradient of the energy in Cartesian coordinates.
       include 'COMMON.MD'
       include 'COMMON.LOCAL'
       include 'COMMON.SPLITELE'
+      include 'COMMON.SCCOR'
       common /srutu/ icall
       dimension ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),
      &  g(maxvar)
@@ -288,10 +291,16 @@ c      call checkintcartgrad
       call geom_to_var(nvar,x)
       if (.not.split_ene) then
         call etotal(energia(0))
+c        do i=1,nres
+c        write (iout,*) "atu?", gloc_sc(1,i,icg),gloc(i,icg)
+c        enddo
         etot=energia(0)
         call enerprint(energia(0))
         call flush(iout)
         write (iout,*) "enter cartgrad"
+c        do i=1,nres
+c        write (iout,*) gloc_sc(1,i,icg)
+c        enddo 
         call flush(iout)
         call cartgrad
         write (iout,*) "exit cartgrad"
@@ -338,6 +347,9 @@ c      call checkintcartgrad
         call zerograd
         call etotal_short(energia(0))
         call enerprint(energia(0))
+c        do i=1,nres
+c        write (iout,*) gloc_sc(1,i,icg)
+c        enddo 
         call flush(iout)
         write (iout,*) "enter cartgrad"
         call flush(iout)
@@ -515,7 +527,7 @@ c-------------------------------------------------------------------------
         be=0.0D0
         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
+        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
index d3e68f7..cd1bd7b 100644 (file)
@@ -1,11 +1,11 @@
 C DO NOT EDIT THIS FILE - IT HAS BEEN GENERATED BY COMPINFO.C
-C 2 5 38
+C 2 5 250
       subroutine cinfo
       include 'COMMON.IOUNITS'
       write(iout,*)'++++ Compile info ++++'
-      write(iout,*)'Version 2.5 build 38'
-      write(iout,*)'compiled Mon Apr  2 08:28:34 2012'
-      write(iout,*)'compiled by adam@matrix.chem.cornell.edu'
+      write(iout,*)'Version 2.5 build 250'
+      write(iout,*)'compiled Thu Jun 14 07:01:36 2012'
+      write(iout,*)'compiled by aks255@matrix.chem.cornell.edu'
       write(iout,*)'OS name:    Linux '
       write(iout,*)'OS release: 2.6.34.9-69.fc13.x86_64 '
       write(iout,*)'OS version:',
@@ -13,7 +13,7 @@ C 2 5 38
       write(iout,*)'flags:'
       write(iout,*)'INSTALL_DIR = /users/software/mpich-1.2.7p1_int...'
       write(iout,*)'FC= ifort'
-      write(iout,*)'OPT =  -O3 -ip -w '
+      write(iout,*)'OPT =  -g -ip -w -CB '
       write(iout,*)'FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include '
       write(iout,*)'FFLAGS1 = -c -w -g -d2 -CA -CB -I$(INSTALL_DIR)...'
       write(iout,*)'FFLAGS2 = -c -w -g -O0 -I$(INSTALL_DIR)/include  '
index 3257273..3979f64 100644 (file)
@@ -483,6 +483,7 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.CONTROL'
       include 'COMMON.TIME1'
       include 'COMMON.MAXGRAD'
+      include 'COMMON.SCCOR'
 #ifdef TIMING
 #ifdef MPI
       time01=MPI_Wtime()
@@ -755,7 +756,6 @@ c      enddo
      &   +wturn3*gel_loc_turn3(i)
      &   +wturn6*gel_loc_turn6(i)
      &   +wel_loc*gel_loc_loc(i)
-     &   +wsccor*gsccor_loc(i)
       enddo
 #ifdef DEBUG
       write (iout,*) "gloc after adding corr"
@@ -5650,7 +5650,7 @@ C Proline-Proline pair is a special case...
      &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
      &  (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
         gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
-c       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
+        write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
       enddo
 ! 6/20/98 - dihedral angle constraints
       edihcnstr=0.0d0
@@ -5765,6 +5765,7 @@ c      do i=1,ndih_constr
         else
           difi=0.0
         endif
+c        write (iout,*) "gloci", gloc(i-3,icg)
 cd        write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii,
 cd     &    rad2deg*phi0(i),  rad2deg*drange(i),
 cd     &    rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
@@ -5835,6 +5836,7 @@ c     lprn=.true.
         enddo
         gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
         gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
+c        write (iout,*) "gloci", gloc(i-3,icg)
       enddo
       return
       end
@@ -5879,18 +5881,21 @@ cc Omicron is flat angle depending on the value of first digit
 c(see comment below)
 
         gloci=0.0D0
-        do intertyp=1,3
+        do intertyp=1,1 !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...SC
-        if (((intertyp.eq.3).and.(itype(i-2).eq.10).or.
-     &      (itype(i-1).eq.10))
-     &    .or. ((intertyp.eq.1).and.(itype(i-2).ne.10))
-     &    .or. ((intertyp.eq.2).and.(itype(i-1).ne.10))) cycle  
-        if ((intertyp.eq.2).and.(i.eq.iphi_start-1)) cycle
-        if ((intertyp.eq.1).and.(i.eq.iphi_end+1)) cycle
+        if (((intertyp.eq.3).and.((itype(i-2).eq.10).or.
+     &      (itype(i-1).eq.10).or.(itype(i-2).eq.21).or.
+     &      (itype(i-1).eq.21)))
+     &    .or. ((intertyp.eq.1).and.((itype(i-2).eq.10)
+     &     .or.(itype(i-2).eq.21)))
+     &    .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or.
+     &      (itype(i-1).eq.21)))) cycle  
+        if ((intertyp.eq.2).and.(i.le.iphi_start-1)) cycle
+        if ((intertyp.eq.1).and.(i.ge.iphi_end+1)) cycle
         do j=1,nterm_sccor(isccori,isccori1)
           v1ij=v1sccor(j,intertyp,isccori,isccori1)
           v2ij=v2sccor(j,intertyp,isccori,isccori1)
@@ -5900,14 +5905,19 @@ c   3 = SC...Ca...Ca...SC
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
         gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
+c        write (iout,*) "WTF",intertyp,i,itype(i),
+c     &gloc_sc(intertyp,i-3,icg)
         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,intertyp,itori,itori1),j=1,6)
      & ,(v2sccor(j,intertyp,itori,itori1),j=1,6)
         gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
+       enddo !intertyp
       enddo
-      enddo
+c        do i=1,nres
+c        write (iout,*) "W@T@F",  gloc_sc(1,i,icg),gloc(i,icg)
+c        enddo
       return
       end
 c----------------------------------------------------------------------------
index 375fcf4..7fec1e8 100644 (file)
@@ -8,6 +8,7 @@
       include 'COMMON.FFIELD'
       include 'COMMON.MD'
       include 'COMMON.IOUNITS'
+      include 'COMMON.SCCOR'
       external ufparm
       integer uiparm(1)
       double precision urparm(1)
@@ -263,10 +264,15 @@ C-------------------------------------------------------------------------
       include 'COMMON.MD'
       include 'COMMON.IOUNITS'
       include 'COMMON.TIME1'
+      include 'COMMON.SCCOR'
 c
 c This subrouting calculates total Cartesian coordinate gradient. 
 c The subroutine chainbuild_cart and energy MUST be called beforehand.
 c
+c        do i=1,nres
+c        write (iout,*) "przed sum_grad", gloc_sc(1,i,icg),gloc(i,icg)
+c        enddo
+
 #ifdef TIMING
       time00=MPI_Wtime()
 #endif
@@ -274,6 +280,9 @@ c
       call sum_gradient
 #ifdef TIMING
 #endif
+c        do i=1,nres
+c        write (iout,*) "checkgrad", gloc_sc(1,i,icg),gloc(i,icg)
+c        enddo     
 cd      write (iout,*) "After sum_gradient"
 cd      do i=1,nres-1
 cd        write (iout,*) i," gradc  ",(gradc(j,i,icg),j=1,3)
@@ -337,6 +346,7 @@ C-------------------------------------------------------------------------
       include 'COMMON.CHAIN'
       include 'COMMON.VAR'
       include 'COMMON.MD'
+      include 'COMMON.SCCOR'
 C
 C Initialize Cartesian-coordinate gradient
 C
@@ -374,6 +384,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 149e86f..b2d378e 100644 (file)
@@ -13,9 +13,15 @@ c-------------------------------------------------------------
       include 'COMMON.INTERACT'
       include 'COMMON.MD'
       include 'COMMON.IOUNITS'
-      
+      include 'COMMON.SCCOR' 
 c   calculating dE/ddc1      
        if (nres.lt.3) goto 18
+c       do i=1,nres
+c c       do intertyp=1,3
+c          write (iout,*) "przed tosyjnymi",i,intertyp,gcart(intertyp,i)
+c     &,gloc_sc(1,i,icg),gloc(i,icg)
+c          enddo
+c       enddo
        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)      
@@ -119,9 +125,13 @@ 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) then
+       if ((itype(1).ne.10).and.(itype(1).ne.21)) then
         do j=1,3
 cc Derviative was calculated for oposite vector of side chain therefore
 c there is "-" sign before gloc_sc
@@ -129,15 +139,16 @@ c there is "-" sign before gloc_sc
      &     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) then
+          if ((itype(2).ne.10).and.(itype(2).ne.21)) 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)
+     &       dtauangle(j,3,2,3)
           endif
        enddo
        endif
-         if (nres.ge.3).and.(itype(3).ne.10) then
+         if ((nres.ge.3).and.(itype(3).ne.10).and.(itype(3).ne.21))
+     & then
          do j=1,3
          gcart(j,1)=gcart(j,1)+gloc_sc(2,1,icg)*dtauangle(j,2,1,4)
          enddo
@@ -148,36 +159,51 @@ 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) then
-           if (itype(1).ne.10) gxcart=(j,2)=gxcart(j,2)+
+         if((itype(2).ne.10).and.(itype(2).ne.21)) 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)) then
+        if ((itype(3).ne.10).and.(nres.ge.3).and.(itype(3).ne.21)) 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(2,1,icg)*dtauangle(j,1,1,4)
+           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(2,1,icg)*dtauangle(j,1,2,4)
+           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) then
+         if ((itype(1).ne.10).and.(itype(1).ne.21)) 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
+         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)
+           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) 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,1i+2)
-          gcart(j,i)=gcart(j,i)+gloc_sc(1,i-1,icg)*
+     &    -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)
@@ -196,6 +222,8 @@ c    If there are more than five residues
           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)*
@@ -207,37 +235,42 @@ c    If there are more than five residues
 c  Setting dE/ddnres-1       
       if(nres.ge.4) then
          do j=1,3
-         if (itype(nres-1).ne.10) then
+         if ((itype(nres-1).ne.10).and.(itype(nres-1).ne.21)) then
          gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(2,nres-3,icg)
      &    *dtauangle(j,2,3,nres)
+          write (iout,*) "gxcart(nres-1)", gloc_sc(2,nres-3,icg),
+     &     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) then
+         if ((itype(nres).ne.10).and.(itype(nres).ne.21)) 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) then
-            gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(2,nres-3,icg)*
-     &   *dtauangle(j,2,3,nres)
+         if ((itype(nres-2).ne.10).and.(itype(nres-2).ne.21)) 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) then
+          if ((itype(nres).ne.10).and.(itype(nres).ne.21)) then
             gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(2,nres-2,icg)*
-          *dtauangle(j,2,2,nres+1)
+     &     dtauangle(j,2,2,nres+1)
+           write (iout,*) "gcart(nres-1)", gloc_sc(2,nres-2,icg),
+     &     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)then
+       if ((nres.ge.3).and.(itype(nres).ne.10))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 8f4a71e..2068031 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),
@@ -51,20 +52,20 @@ c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
 #else
       do i=3,nres
 #endif
-      if (itype(i-1).ne.10) then
+      if ((itype(i-1).ne.10).and.(itype(i-1).ne.21)) 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)+
+          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+nres)
+     &    +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
@@ -73,8 +74,9 @@ CC Looks messy but better than if in loop
      &    vbld(i)
           domicron(j,2,1,i)=-1/sint2*dcosomicron(j,2,1,i)
           dcosomicron(j,2,2,i)=-(dc_norm(j,i-1)
-     &     +cost1*(-dc_norm(j,i-1+nres)))/
-     &    vbld(i+nres)
+     &     +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
@@ -153,6 +155,7 @@ Calculate derivative of Tauangle
 #else
       do i=3,nres
 #endif
+       if ((itype(i-2).eq.21).or.(itype(i-2).eq.10)) cycle
 cc dtauangle(j,intertyp,dervityp,residue number)
 cc INTERTYP=1 SC...Ca...Ca..Ca
 c the conventional case
@@ -164,6 +167,7 @@ c the conventional case
         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)
@@ -171,13 +175,11 @@ c the conventional case
         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
-         do j=1,3
-         dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
-         enddo
          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)
@@ -185,37 +187,42 @@ c    Obtaining the gamma derivatives from sine derivative
             ctgt=cost/sint
             ctgt1=cost1/sint1
             cosg_inv=1.0d0/cosg
-            dsintau(j,1,1,i)=-sing*ctgt1*domicron(j,2,1,i-1)
-     &-(fac0*vp1(j)+sing*(-dc_norm(j,i-2+nres)))
+            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,2,i-1)+ctgt*dtheta(j,1,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*dsinphi(j,1,3,i)
+            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,1,i-1)+fac3*
-     &     dcosomicron(j,2,1,i-1)-fac0*(dc_norm(j,i-1)-scalp*
-     &     (-dc_norm(j,i-2+nres)))/vbld(i-2+nres)
+           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,2,i-1)+fac2*
-     &     dcostheta(j,1,i)+fac3*dcosomicron(j,2,2,i-1)+fac4*
+           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                                                                                            
+        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
@@ -223,11 +230,12 @@ CC Second case Ca...Ca...Ca...SC
 #else
       do i=4,nres
 #endif
+       if ((itype(i-1).eq.21).or.(itype(i-1).eq.10)) 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)
+        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
@@ -251,17 +259,22 @@ c    Obtaining the gamma derivatives from sine derivative
             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)
+     &        +(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)
-            dphi(j,2,i)=cosg_inv*dsintau(j,2,2,i)
+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)
+     &       +(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*dsinphi(j,2,3,i)
+            dtauangle(j,2,3,i)=cosg_inv*dsintau(j,2,3,i)
          enddo                                          
 c   Obtaining the gamma derivatives from cosine derivative
         else
@@ -273,11 +286,12 @@ c   Obtaining the gamma derivatives from cosine derivative
            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)
-           dtauanlge(j,2,2,i)=-1/sing*dcostau(j,2,2,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*
-     &     dcostheta(j,1,2,i)-fac0*(dc_norm(j,i-3)-scalp*
+     &     dcosomicron(j,1,2,i)-fac0*(dc_norm(j,i-3)-scalp*
      &     dc_norm(j,i-1+nres))/vbld(i-1+nres)
-           dtauanlge(j,2,3,i)=-1/sing*dcosphi(j,3,i)
+           dtauangle(j,2,3,i)=-1/sing*dcostau(j,2,3,i)
+        write(iout,*) i,j,"else", dtauangle(j,2,3,i) 
          enddo
         endif                                                                                            
       enddo
@@ -291,6 +305,8 @@ CCC third case SC...Ca...Ca...SC
       do i=3,nres
 #endif
 c the conventional case
+      if ((itype(i-1).eq.21).or.(itype(i-1).eq.10).or.
+     &(itype(i-2).eq.21).or.(itype(i-2).eq.10)) cycle
         sint=dsin(omicron(1,i))
         sint1=dsin(omicron(2,i-1))
         sing=dsin(tauangle(3,i))
@@ -331,7 +347,7 @@ c Bug fixed 3/24/05 (AL)
      &        +(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)
-            dphi(j,3,3,i)=cosg_inv*dsintau(j,3,3,i)
+            dtauangle(j,3,3,i)=cosg_inv*dsintau(j,3,3,i)
          enddo                                          
 c   Obtaining the gamma derivatives from cosine derivative
         else
index be1a193..1968235 100644 (file)
@@ -562,15 +562,17 @@ C
 c      write (iout,*) 'ntortyp',ntortyp
       maxinter=3
 cc maxinter is maximum interaction sites
-      do l=1,maxinter
+      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,i,j),v2sccor(k,i,j) 
-            v0ijsccor=v0ijsccor+si*v1sccor(k,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)
@@ -592,7 +594,7 @@ cc maxinter is maximum interaction sites
             write (iout,*) 'ityp',i,' jtyp',j
             write (iout,*) 'Fourier constants'
             do k=1,nterm_sccor(i,j)
-             write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,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)