debug w poszukiwaniu bledu w ostatniej reszcie
authorAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Fri, 19 Dec 2014 12:29:56 +0000 (13:29 +0100)
committerAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Fri, 19 Dec 2014 12:29:56 +0000 (13:29 +0100)
source/unres/src_MD-M/COMMON.TORSION
source/unres/src_MD-M/energy_p_new_barrier.F
source/unres/src_MD-M/gradient_p.F
source/unres/src_MD-M/int_to_cart.f
source/unres/src_MD-M/parmread.F
source/unres/src_MD-M/stochfric.F

index 2a887a6..e0f9c5d 100644 (file)
@@ -10,7 +10,7 @@ C Torsional constants of the rotation about virtual-bond dihedral angles
      &    nterm(maxtor,maxtor),
      &    nlor(maxtor,maxtor), 
      &    nterm_old,
-     &    itortyp(ntyp),
+     &    itortyp(ntyp1),
      &    ntortyp
 C 6/23/01 - constants for double torsionals
       double precision v1c,v1s,v2c,v2s
index b9d17fe..0c2f5ed 100644 (file)
@@ -85,7 +85,7 @@ C FG slaves receive the WEIGHTS array
         time_Bcastw=time_Bcastw+MPI_Wtime()-time00
 c        call chainbuild_cart
       endif
-c      print *,'Processor',myrank,' calling etotal ipot=',ipot
+C      print *,'Processor',myrank,' calling etotal ipot=',ipot
 c      print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
 #else
 c      if (modecalc.eq.12.or.modecalc.eq.14) then
@@ -98,6 +98,7 @@ c      endif
 C 
 C Compute the side-chain and electrostatic interaction energy
 C
+C      write(iout,*) "zaczynam liczyc energie"
       goto (101,102,103,104,105,106) ipot
 C Lennard-Jones potential.
   101 call elj(evdw)
@@ -117,10 +118,13 @@ C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
       goto 107
 C Soft-sphere potential
   106 call e_softsphere(evdw)
+C      write(iout,*) "skonczylem ipoty"
+
 C
 C Calculate electrostatic (H-bonding) energy of the main chain.
 C
   107 continue
+C           write(iout,*) "skonczylem ipoty"
 cmc
 cmc Sep-06: egb takes care of dynamic ss bonds too
 cmc
@@ -442,7 +446,6 @@ cMS$ATTRIBUTES C ::  proc_proc
 #endif
       double precision gradbufc(3,maxres),gradbufx(3,maxres),
      &  glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
-#endif
       include 'COMMON.SETUP'
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
@@ -2338,6 +2341,7 @@ C
         endif
 c        if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
         if (i.gt. nnt+2 .and. i.lt.nct+2) then
+c        write(iout,*) (itype(i-2))
           iti = itortyp(itype(i-2))
         else
           iti=ntortyp+1
@@ -2718,7 +2722,7 @@ C
       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
-     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
+     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),eel_loc_ij
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
      &    num_conti,j1,j2
@@ -2761,6 +2765,7 @@ c        call vec_and_deriv
         time01=MPI_Wtime()
 #endif
         call set_matrices
+c        write (iout,*) "after set matrices"
 #ifdef TIMING
         time_mat=time_mat+MPI_Wtime()-time01
 #endif
@@ -2797,6 +2802,7 @@ c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
 C
 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
 C
+c      write(iout,*) "przed turnem3 loop"
       do i=iturn3_start,iturn3_end
         if (itype(i).eq.21 .or. itype(i+1).eq.21 
      &  .or. itype(i+2).eq.21 .or. itype(i+3).eq.21) cycle
@@ -2889,7 +2895,7 @@ C-------------------------------------------------------------------------------
       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
-     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
+     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),a22,a23,a32,a33
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
      &    num_conti,j1,j2
@@ -3659,7 +3665,7 @@ c        write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2
         iti1=itortyp(itype(i+1))
         iti2=itortyp(itype(i+2))
         iti3=itortyp(itype(i+3))
-c        write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3
+        write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3
         call transpose2(EUg(1,1,i+1),e1t(1,1))
         call transpose2(Eug(1,1,i+2),e2t(1,1))
         call transpose2(Eug(1,1,i+3),e3t(1,1))
@@ -5646,7 +5652,7 @@ C Set lprn=.true. for debugging
       lprn=.false.
 c     lprn=.true.
       etors_d=0.0D0
-c      write(iout,*) "a tu??"
+C      write(iout,*) "a tu??"
       do i=iphid_start,iphid_end
         if (itype(i-2).eq.21 .or. itype(i-1).eq.21
      &      .or. itype(i).eq.21 .or. itype(i+1).eq.21) cycle
index 2c670f2..b5ba82a 100644 (file)
@@ -292,9 +292,11 @@ c If performing constraint dynamics, add the gradients of the constraint energy
          do i=1,nres-3
            gloc(i,icg)=gloc(i,icg)+dugamma(i)
          enddo
+      write(iout,*) "TU JESTEM?"
          do i=1,nres-2
            gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
          enddo
+        write(iout,*) "TU JESTEM?"
       endif 
 #ifdef TIMING
       time01=MPI_Wtime()
index 24ce0bc..e178afe 100644 (file)
@@ -24,15 +24,16 @@ c   calculating dE/ddc1
      &    gloc(ialph(2,1)+nside,icg)*domega(j,1,2)             
         endif
        enddo
+C      write (iout,*) "????A CO TO??"
 c     Calculating the remainder of dE/ddc2
        do j=1,3
          gcart(j,2)=gcart(j,2)+gloc(1,icg)*dphi(j,2,4)+
      &  gloc(nres-2,icg)*dtheta(j,2,3)+gloc(nres-1,icg)*dtheta(j,1,4)
-        if(itype(2).ne.10) then
+        if((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then
           gcart(j,2)=gcart(j,2)+gloc(ialph(2,1),icg)*dalpha(j,2,2)+
      &    gloc(ialph(2,1)+nside,icg)*domega(j,2,2)
         endif
-               if(itype(3).ne.10) then
+               if((itype(3).ne.10).and.(itype(3).ne.21)) then
          gcart(j,2)=gcart(j,2)+gloc(ialph(3,1),icg)*dalpha(j,1,3)+
      &    gloc(ialph(3,1)+nside,icg)*domega(j,1,3)
         endif
@@ -56,6 +57,7 @@ c  If there are only five residues
          endif
        enddo
        endif
+C      write (iout,*) "Poniezej blad??",ialph(3,1),nside,ialph(4,1)
 c    If there are more than five residues
       if(nres.gt.5) then                          
         do i=3,nres-3
@@ -64,30 +66,31 @@ c    If there are more than five residues
      &    +gloc(i-1,icg)*dphi(j,2,i+2)+
      &    gloc(i,icg)*dphi(j,1,i+3)+gloc(nres+i-4,icg)*dtheta(j,2,i+1)+
      &    gloc(nres+i-3,icg)*dtheta(j,1,i+2)
-          if(itype(i).ne.10) then
+          if((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
            gcart(j,i)=gcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,2,i)+
      &     gloc(ialph(i,1)+nside,icg)*domega(j,2,i)
           endif
-          if(itype(i+1).ne.10) then
+          if((itype(i+1).ne.10).and.(itype(i+1).ne.ntyp1)) then
            gcart(j,i)=gcart(j,i)+gloc(ialph(i+1,1),icg)*dalpha(j,1,i+1)
      &     +gloc(ialph(i+1,1)+nside,icg)*domega(j,1,i+1)
           endif
          enddo
         enddo
       endif    
-c  Setting dE/ddnres-2       
+c  Setting dE/ddnres-2      
+      write(iout,*) "ATUCHUJ?" 
       if(nres.gt.5) then
          do j=1,3
            gcart(j,nres-2)=gcart(j,nres-2)+gloc(nres-4,icg)*
      &    dphi(j,3,nres-1)+gloc(nres-3,icg)*dphi(j,2,nres)
      &     +gloc(2*nres-6,icg)*
      &     dtheta(j,2,nres-1)+gloc(2*nres-5,icg)*dtheta(j,1,nres)
-          if(itype(nres-2).ne.10) then
+          if((itype(nres-2).ne.10).and.(itype(nres-2).ne.ntyp1)) then
               gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-2,1),icg)*
      &       dalpha(j,2,nres-2)+gloc(ialph(nres-2,1)+nside,icg)*
      &        domega(j,2,nres-2)
           endif
-          if(itype(nres-1).ne.10) then
+          if((itype(nres-1).ne.10).and.(itype(nres-1).ne.ntyp1)) then
              gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-1,1),icg)*
      &      dalpha(j,1,nres-1)+gloc(ialph(nres-1,1)+nside,icg)*
      &       domega(j,1,nres-1)
@@ -98,7 +101,7 @@ c  Settind dE/ddnres-1
        do j=1,3
         gcart(j,nres-1)=gcart(j,nres-1)+gloc(nres-3,icg)*dphi(j,3,nres)+
      & gloc(2*nres-5,icg)*dtheta(j,2,nres)
-        if(itype(nres-1).ne.10) then
+        if((itype(nres-1).ne.10).and.(itype(nres-1).ne.ntyp1)) then
           gcart(j,nres-1)=gcart(j,nres-1)+gloc(ialph(nres-1,1),icg)*
      &   dalpha(j,2,nres-1)+gloc(ialph(nres-1,1)+nside,icg)*
      &    domega(j,2,nres-1)
@@ -106,13 +109,15 @@ c  Settind dE/ddnres-1
         enddo
 c   The side-chain vector derivatives
         do i=2,nres-1
-         if(itype(i).ne.10 .and. itype(i).ne.21) then  
+         if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then       
             do j=1,3   
               gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i)
      &        +gloc(ialph(i,1)+nside,icg)*domega(j,3,i)
             enddo
          endif     
-       enddo                                                                                                                                                   
+       enddo
+       write(iout,*) "TU DOCHODZE"
+                                                                                                                                                       
 c----------------------------------------------------------------------
 C INTERTYP=1 SC...Ca...Ca...Ca
 C INTERTYP=2 Ca...Ca...Ca...SC
index 4dd6b0a..9bc2d20 100644 (file)
@@ -397,6 +397,7 @@ C Read torsional parameters
 C
       read (itorp,*,end=113,err=113) ntortyp
       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
+      itortyp(ntyp1)=ntortyp+1
 c      write (iout,*) 'ntortyp',ntortyp
       do i=1,ntortyp
        do j=1,ntortyp
@@ -609,7 +610,7 @@ cc maxinter is maximum interaction sites
             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
      &(1+vlor3sccor(k,i,j)**2)
           enddo
-          v0sccor(i,j)=v0ijsccor
+          v0sccor(l,i,j)=v0ijsccor
         enddo
       enddo
       enddo
index 8d1cf2c..76aa40a 100644 (file)
         ind=ind+3
       enddo
       do i=nnt,nct
-<<<<<<< HEAD
         if (itype(i).ne.10 .and. itype(i).ne.21) then
-=======
-        if (itype(i).ne.10) then
->>>>>>> 2acc991... Wprowadznie potencjalow SC-COR do multichain oraz ich pseudosymetrii dla
           do j=1,3
             d_t_work(ind+j)=d_t(j,i+nres)
           enddo
         ind=ind+3
       enddo
       do i=nnt,nct
-<<<<<<< HEAD
         if (itype(i).ne.10 .and. itype(i).ne.21) then
-=======
-        if (itype(i).ne.10) then
->>>>>>> 2acc991... Wprowadznie potencjalow SC-COR do multichain oraz ich pseudosymetrii dla
           do j=1,3
             friction(j,i+nres)=fric_work(ind+j)
           enddo
@@ -242,11 +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
-<<<<<<< HEAD
         if (itype(i).ne.10 .and. itype(i).ne.21) then
-=======
-        if (itype(i).ne.10) then
->>>>>>> 2acc991... Wprowadznie potencjalow SC-COR do multichain oraz ich pseudosymetrii dla
           do j=1,3
             stochforc(j,i+nres)=force(j,i+nres)
           enddo
@@ -264,11 +252,7 @@ c Compute the stochastic forces acting on virtual-bond vectors.
         ind=ind+3
       enddo
       do i=nnt,nct
-<<<<<<< HEAD
         if (itype(i).ne.10 .and. itype(i).ne.21) then
-=======
-        if (itype(i).ne.10) then
->>>>>>> 2acc991... Wprowadznie potencjalow SC-COR do multichain oraz ich pseudosymetrii dla
           do j=1,3
             stochforcvec(ind+j)=stochforc(j,i+nres)
           enddo