Wprowadznie potencjalow SC-COR do multichain oraz ich pseudosymetrii dla
authorAdam Kazimierz Sieradzan <adasko@sun1.chem.univ.gda.pl>
Thu, 4 Oct 2012 14:52:15 +0000 (10:52 -0400)
committerAdam Kazimierz Sieradzan <adasko@sun1.chem.univ.gda.pl>
Thu, 4 Oct 2012 14:52:15 +0000 (10:52 -0400)
D aminokwasow

21 files changed:
bin/unres/MD-M/unres_Tc_procor_oldparm_em64-D-symetr.exe
source/unres/src_MD-M/COMMON.CHAIN
source/unres/src_MD-M/COMMON.LOCAL
source/unres/src_MD-M/COMMON.SCCOR
source/unres/src_MD-M/COMMON.VAR
source/unres/src_MD-M/DIMENSIONS
source/unres/src_MD-M/MD.F
source/unres/src_MD-M/Makefile
source/unres/src_MD-M/checkder_p.F
source/unres/src_MD-M/cinfo.f
source/unres/src_MD-M/energy_p_new_barrier.F
source/unres/src_MD-M/gradient_p.F
source/unres/src_MD-M/initialize_p.F
source/unres/src_MD-M/int_to_cart.f
source/unres/src_MD-M/intcartderiv.F
source/unres/src_MD-M/kinetic_lesyng.f
source/unres/src_MD-M/lagrangian_lesyng.F
source/unres/src_MD-M/moments.f
source/unres/src_MD-M/parmread.F
source/unres/src_MD-M/stochfric.F
source/unres/src_MD/COMMON.SCCOR

index 0cc9b64..71157fe 100755 (executable)
Binary files a/bin/unres/MD-M/unres_Tc_procor_oldparm_em64-D-symetr.exe and b/bin/unres/MD-M/unres_Tc_procor_oldparm_em64-D-symetr.exe differ
index ea1c17a..71db8f7 100644 (file)
@@ -1,9 +1,10 @@
       integer nres,nsup,nstart_sup,nz_start,nz_end,iz_sc,
      &  nres0,nstart_seq,chain_length,iprzes,tabperm,nperm
       double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm,t,r,
-     & prod,rt,dc_work,cref,crefjlee,chain_rep
+     & prod,rt,dc_work,cref,crefjlee,chain_rep,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 e1e13f4..c8d686f 100644 (file)
@@ -34,7 +34,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),
@@ -48,7 +48,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 5217de7..e29cb4c 100644 (file)
@@ -1,6 +1,18 @@
-C Parameters of the SCCOR term
-      double precision v1sccor,v2sccor
-      integer nterm_sccor
-      common/torsion/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 d9992af..5131f4e 100644 (file)
@@ -59,7 +59,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 a28d660..f65a527 100644 (file)
@@ -1628,7 +1628,7 @@ c if the friction coefficients do not depend on surface area
           stdforcp(i)=stdfp*dsqrt(gamp)
         enddo
         do i=nnt,nct
-          stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
+          stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(iabs(itype(i))))
         enddo
       endif
 c Open the pdb file for snapshotshots
index e09600e..35006b3 100644 (file)
@@ -1,7 +1,7 @@
 CPPFLAGS = -DLINUX -DUNRES -DMP -DMPI \
            -DPGI -DSPLITELE -DISNAN -DAMD64 \
            -DPROCOR -DLANG0 \
-       -DCRYST_BOND -DCRYST_THETA -DCRYST_SC 
+       -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DSCCORPDB 
 ## -DPROCOR
 ## -DMOMENT
 #-DCO_BIAS
@@ -14,8 +14,9 @@ INSTALL_DIR = /users/software/mpich-1.2.7p1_intel-10.1_em64_ssh
 #FC= /usr/local/opt/intel/compiler60/ia32/bin/ifc
 FC= ifort
 
-OPT =  -O3 -ip -w
-#OPT = -g -CB
+#OPT =  -O3 -ip -w
+OPT = -g -CB
+#OPT = -g
 CFLAGS = -DSGI -c
 
 FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include
index 26854e6..0539e48 100644 (file)
@@ -513,7 +513,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 f46c7f9..4e7d75d 100644 (file)
@@ -1,10 +1,10 @@
 C DO NOT EDIT THIS FILE - IT HAS BEEN GENERATED BY COMPINFO.C
-C 2 3 3404
+C 2 3 3438
       subroutine cinfo
       include 'COMMON.IOUNITS'
       write(iout,*)'++++ Compile info ++++'
-      write(iout,*)'Version 2.3 build 3404'
-      write(iout,*)'compiled Thu Aug 16 11:51:48 2012'
+      write(iout,*)'Version 2.3 build 3438'
+      write(iout,*)'compiled Thu Oct  4 09:08:03 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 '
@@ -14,10 +14,10 @@ C 2 3 3404
       write(iout,*)'CPPFLAGS = -DLINUX -DUNRES -DMP -DMPI \\'
       write(iout,*)'           -DPGI -DSPLITELE -DISNAN -DAMD64 \\'
       write(iout,*)'           -DPROCOR -DLANG0 \\'
-      write(iout,*)'   -DCRYST_BOND -DCRYST_THETA -DCRYST_SC '
+      write(iout,*)'   -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DSCCORP...'
       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 -CB'
       write(iout,*)'CFLAGS = -DSGI -c'
       write(iout,*)'FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include'
       write(iout,*)'FFLAGS1 = -c -w -g -d2 -CA -CB -I$(INSTALL_DIR)...'
index 0ce1a7b..38c5306 100644 (file)
@@ -434,7 +434,7 @@ cMS$ATTRIBUTES C ::  proc_proc
 #ifdef MPI
       include 'mpif.h'
       double precision gradbufc(3,maxres),gradbufx(3,maxres),
-     &  glocbuf(4*maxres),gradbufc_sum(3,maxres)
+     &  glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
 #endif
       include 'COMMON.SETUP'
       include 'COMMON.IOUNITS'
@@ -447,6 +447,7 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.CONTROL'
       include 'COMMON.TIME1'
       include 'COMMON.MAXGRAD'
+      include 'COMMON.SCCOR'
 #ifdef TIMING
       time01=MPI_Wtime()
 #endif
@@ -689,7 +690,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"
@@ -708,6 +708,21 @@ c      enddo
         do i=1,4*nres
           glocbuf(i)=gloc(i,icg)
         enddo
+#define DEBUG
+#ifdef DEBUG
+      write (iout,*) "gloc_sc before reduce"
+      do i=1,nres
+       do j=1,1
+        write (iout,*) i,j,gloc_sc(j,i,icg)
+       enddo
+      enddo
+#endif
+#undef DEBUG
+        do i=1,nres
+         do j=1,3
+          gloc_scbuf(j,i)=gloc_sc(j,i,icg)
+         enddo
+        enddo
         time00=MPI_Wtime()
         call MPI_Barrier(FG_COMM,IERR)
         time_barrier_g=time_barrier_g+MPI_Wtime()-time00
@@ -719,6 +734,19 @@ c      enddo
         call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres,
      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
         time_reduce=time_reduce+MPI_Wtime()-time00
+        call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres,
+     &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
+        time_reduce=time_reduce+MPI_Wtime()-time00
+#define DEBUG
+#ifdef DEBUG
+      write (iout,*) "gloc_sc after reduce"
+      do i=1,nres
+       do j=1,1
+        write (iout,*) i,j,gloc_sc(j,i,icg)
+       enddo
+      enddo
+#endif
+#undef DEBUG
 #ifdef DEBUG
       write (iout,*) "gloc after reduce"
       do i=1,4*nres
@@ -5706,29 +5734,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
-        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1) cycle
+      do i=itau_start,itau_end
         esccor_ii=0.0D0
-        itori=iabs(itype(i-2))
-        itori1=iabs(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 effd955..2c670f2 100644 (file)
@@ -345,6 +345,7 @@ C-------------------------------------------------------------------------
       include 'COMMON.CHAIN'
       include 'COMMON.VAR'
       include 'COMMON.MD'
+      include 'COMMON.SCCOR'
 C
 C Initialize Cartesian-coordinate gradient
 C
@@ -380,6 +381,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 aa5dfb9..ed4fb21 100644 (file)
@@ -591,6 +591,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
@@ -1119,6 +1122,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
index 9208a80..f413622 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 33d4a50..369a4f0 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),
@@ -53,7 +54,42 @@ c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
           if (itype(i-1).ne.ntyp1) dtheta(j,2,i)=-dcostheta(j,2,i)/sint
         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
@@ -123,6 +159,228 @@ 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 44ea85b..db959b3 100644 (file)
@@ -44,7 +44,7 @@ c to the velocities of the first Calpha.
         incr(j)=d_t(j,0)
       enddo
       do i=nnt,nct
-        iti=itype(i)
+        iti=iabs(itype(i))
         if (itype(i).eq.10) then
           do j=1,3
             v(j)=incr(j)
index 9611f90..f8834ea 100644 (file)
@@ -217,12 +217,12 @@ c  Diagonal elements of the dX part of A and the respective friction coefficient
         ind=ind+1
         ii = ind+m
         iti=itype(i)
-        massvec(ii)=msc(iti)
+        massvec(ii)=msc(iabs(iti))
         if (iti.ne.10 .and. iti.ne.ntyp1) then
           ind1=ind1+1
           ii1= ind1+m1
           A(ii,ii1)=1.0d0
-          Gmat(ii1,ii1)=ISC(iti)
+          Gmat(ii1,ii1)=ISC(iabs(iti))
         endif
       enddo
 c  Off-diagonal elements of the dX part of A
index 5992773..983ce36 100644 (file)
@@ -43,10 +43,10 @@ c   calculating the center of the mass of the protein
         M_SC=0.0d0                             
         do i=nnt,nct
            iti=iabs(itype(i))           
-          M_SC=M_SC+msc(iti)
+          M_SC=M_SC+msc(iabs(iti))
            inres=i+nres
            do j=1,3
-            cm(j)=cm(j)+msc(iti)*c(j,inres)        
+            cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)          
            enddo
         enddo
         do j=1,3
@@ -71,12 +71,12 @@ c   calculating the center of the mass of the protein
            do j=1,3
              pr(j)=c(j,inres)-cm(j)        
            enddo
-          Im(1,1)=Im(1,1)+msc(iti)*(pr(2)*pr(2)+pr(3)*pr(3))
-          Im(1,2)=Im(1,2)-msc(iti)*pr(1)*pr(2)
-          Im(1,3)=Im(1,3)-msc(iti)*pr(1)*pr(3)
-          Im(2,3)=Im(2,3)-msc(iti)*pr(2)*pr(3) 
-          Im(2,2)=Im(2,2)+msc(iti)*(pr(3)*pr(3)+pr(1)*pr(1))
-          Im(3,3)=Im(3,3)+msc(iti)*(pr(1)*pr(1)+pr(2)*pr(2))              
+          Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3))
+          Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2)
+          Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3)
+          Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3)   
+          Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1))
+          Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2))                
         enddo
           
         do i=nnt,nct-1
@@ -262,7 +262,7 @@ c  Calculate the angular momentum
 c         write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
 c     &     " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
          do j=1,3
-            L(j)=L(j)+msc(iti)*vp(j)
+            L(j)=L(j)+msc(iabs(iti))*vp(j)
          enddo
 c         write (iout,*) "L",(l(j),j=1,3)
          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
index 9e1614b..33f017c 100644 (file)
@@ -600,32 +600,129 @@ C Martix of D parameters for two dimesional fourier series
       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,ntyp
-       do j=1,ntyp
-          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)
+      do i=-ntyp,-1
+        isccortyp(i)=-isccortyp(-i)
+      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
+          nterm_sccor(-i,j)=nterm_sccor(i,j)
+          nterm_sccor(-i,-j)=nterm_sccor(i,j)
+          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)
+            if (j.eq.iscprol) then
+             if (i.eq.isccortyp(10)) then
+             v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
+             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+             else
+             v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
+     &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
+             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
+     &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
+             v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
+             v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
+             v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
+             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
+             endif
+            else
+             if (i.eq.isccortyp(10)) then
+             v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
+             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+             else
+               if (j.eq.isccortyp(10)) then
+             v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
+             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
+               else
+             v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
+             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+             v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
+             v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
+             v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
+             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
+                endif
+               endif
+            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,ntyp
-         do j=1,ntyp
+        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)
+            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
       endif
+
 C
 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 C         interaction energy of the Gly, Ala, and Pro prototypes.
index 6427551..cda93c3 100644 (file)
@@ -39,7 +39,7 @@
         ind=ind+3
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+        if (itype(i).ne.10) then
           do j=1,3
             d_t_work(ind+j)=d_t(j,i+nres)
           enddo
@@ -68,7 +68,7 @@
         ind=ind+3
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+        if (itype(i).ne.10) then
           do j=1,3
             friction(j,i+nres)=fric_work(ind+j)
           enddo
@@ -185,7 +185,11 @@ c-----------------------------------------------------
       enddo
       x=0.0d0  
 
+#ifdef MPI
       time00=MPI_Wtime()
+#else
+      time00=tcpu()
+#endif
 c Compute the stochastic forces acting on bodies. Store in force.
       do i=nnt,nct-1
         sig=stdforcp(i)
@@ -203,7 +207,11 @@ c Compute the stochastic forces acting on bodies. Store in force.
           force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
         enddo
       enddo
+#ifdef MPI
       time_fsample=time_fsample+MPI_Wtime()-time00
+#else
+      time_fsample=time_fsample+tcpu()-time00
+#endif
 c Compute the stochastic forces acting on virtual-bond vectors.
       do j=1,3
         ff(j)=0.0d0
@@ -226,7 +234,7 @@ c Compute the stochastic forces acting on virtual-bond vectors.
         stochforc(j,0)=ff(j)+force(j,nnt+nres)
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+        if (itype(i).ne.10) then
           do j=1,3
             stochforc(j,i+nres)=force(j,i+nres)
           enddo
@@ -244,7 +252,7 @@ c Compute the stochastic forces acting on virtual-bond vectors.
         ind=ind+3
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+        if (itype(i).ne.10) then
           do j=1,3
             stochforcvec(ind+j)=stochforc(j,i+nres)
           enddo
@@ -309,7 +317,9 @@ c Compute the stochastic forces acting on virtual-bond vectors.
 c------------------------------------------------------------------
       subroutine setup_fricmat
       implicit real*8 (a-h,o-z)
+#ifdef MPI
       include 'mpif.h'
+#endif
       include 'DIMENSIONS'
       include 'COMMON.VAR'
       include 'COMMON.CHAIN'
@@ -355,7 +365,6 @@ c  Load the friction coefficients corresponding to peptide groups
 c  Load the friction coefficients corresponding to side chains
       m=nct-nnt
       ind=0
-      gamsc(ntyp1)=1.0d0
       do i=nnt,nct
         ind=ind+1
         ii = ind+m
@@ -478,22 +487,39 @@ c      enddo
       if (nfgtasks.gt.1) then
         if (fg_rank.eq.0) then
 c The matching BROADCAST for fg processors is called in ERGASTULUM
+#ifdef MPI
           time00=MPI_Wtime()
+#else
+          time00=tcpu()
+#endif
           call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
+#ifdef MPI
           time_Bcast=time_Bcast+MPI_Wtime()-time00
+#else
+          time_Bcast=time_Bcast+tcpu()-time00
+#endif
 c          print *,"Processor",myrank,
 c     &       " BROADCAST iorder in SETUP_FRICMAT"
         endif
 c      licznik=licznik+1
 c        write (iout,*) "setup_fricmat licznik",licznik
+#ifdef MPI
         time00=MPI_Wtime()
+#else
+        time00=tcpu()
+#endif
 c Scatter the friction matrix
         call MPI_Scatterv(fricmat(1,1),nginv_counts(0),
      &    nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),
      &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
-        time_scatter=time_scatter+MPI_Wtime()-time00
 #ifdef TIMING
+#ifdef MPI
+        time_scatter=time_scatter+MPI_Wtime()-time00
         time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
+#else
+        time_scatter=time_scatter+tcpu()-time00
+        time_scatter_fmat=time_scatter_fmat+tcpu()-time00
+#endif
 #endif
         do i=1,dimen
           do j=1,2*my_ng_count
@@ -532,7 +558,7 @@ c
       include 'COMMON.NAMES'
       double precision radius(maxres2),gamvec(maxres6)
       parameter (twosix=1.122462048309372981d0)
-      logical lprn /.true./
+      logical lprn /.false./
 c
 c     determine new friction coefficients every few SD steps
 c
index 6e39f93..e29cb4c 100644 (file)
@@ -4,9 +4,9 @@ cc Parameters of the SCCOR term
      &                 dcostau,dsintau,dtauangle,dcosomicron,
      &                 domicron
       integer nterm_sccor,isccortyp,nsccortyp,nlor_sccor
-      common/sccor/v1sccor(maxterm_sccor,3,20,20),
-     &    v2sccor(maxterm_sccor,3,20,20),
-     &    v0sccor(maxterm_sccor,20),
+      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),