in debug mode on

[unres4.git] / source / unres / REMD.f90

diff --git a/source/unres/REMD.f90 b/source/unres/REMD.f90
index 1e294ad..358cf78 100644 (file)
--- a/source/unres/REMD.f90
+++ b/source/unres/REMD.f90
@@ -47,7 +47,7 @@
 #ifdef CHECK5DSOL
        real(kind=8) :: rscheck(dimen),rsold(dimen)
 #endif
-       logical :: lprn = .true.
+       logical :: lprn = .false.
 !el       common /cipiszcze/ itime
        itime = itt_comm
 
@@ -62,14 +62,16 @@
         enddo
         write (iout,*) "Potential forces sidechain"
         do i=nnt,nct
-          if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) &
+          mnum=molnum(i)
+          if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.eq.5)&
             write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3)
         enddo
       endif
        do j=1,3
          ind=1
          do i=nnt,nct
-           if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1)then
+           mnum=molnum(i)     
+           if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1 .or. mnum.eq.5)then
            rs(ind)=-gcart(j,i)-gxcart(j,i)
             ind=ind+1
           else
@@ -192,7 +194,8 @@
       if (lprn) write (iout,*) "Potential forces sidechain"
       do i=nnt,nct
           mnum=molnum(i)
-        if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
+        if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+          .and.mnum.ne.5) then
           if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
              i,(-gxcart(j,i),j=1,3)
           do j=1,3
@@ -217,7 +220,8 @@
       do i=nnt,nct
        mnum=molnum(i)
 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
-         if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
+         if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+          .and.mnum.ne.5) then
           do j=1,3
             ind=ind+1
             d_a(j,i+nres)=d_a_work(ind)
@@ -296,7 +300,7 @@
 !#endif
 !      include 'COMMON.IOUNITS'
 !      include 'COMMON.TIME1'
-      logical :: lprn = .true.
+      logical :: lprn = .false.
       logical :: osob
 #ifndef FIVEDIAG
       real(kind=8),dimension(2*nres,2*nres) :: Bmat,matmult
@@ -311,16 +315,18 @@
       real(kind=8),dimension(8*6*nres) :: work !(8*maxres6)
       integer,dimension(6*nres) :: iwork       !(maxres6) maxres6=6*maxres
 !el      common /przechowalnia/ Gcopy,Ghalf
-      real(kind=8) :: coeff
+      real(kind=8) :: coeff,mscab
       integer :: i,j,ind,ind1,k,l,ii,jj,m,m1,ii1,iti,nres2,ierr,nind,mark
       real(kind=8) :: ip4
       integer :: iz,mnum
+      print *,"just entered"
       relfeh=1.0d-14
       nres2=2*nres
-
+      print *,"FIVEDIAG"
       write (iout,*) "before FIVEDIAG"
 #ifndef FIVEDIAG
       write (iout,*) "ALLOCATE"
+      print *,"ALLOCATE"
       if(.not.allocated(Gcopy)) allocate(Gcopy(nres2,nres2)) !(maxres2,maxres2)
       if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
 #endif
@@ -374,7 +380,9 @@
        
         ind=1
         do i=nnt,nct
-       if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1) then
+        mnum=molnum(i)
+       if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1 &
+          .and.mnum.eq.5) then
           DU1(ind)=-isc(iabs(itype(i,1)))
           DU1(ind+1)=0.0d0
          ind=ind+2
@@ -513,6 +521,8 @@
         ind=ind+1
         ind1=ind1+1
         coeff=0.25d0*IP(mnum)
+        if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+        print *,i,coeff,mp(mnum)
         massvec(ind1)=mp(mnum)
         Gmat(ind,ind)=coeff
         print *,"i",mp(mnum)
@@ -539,7 +549,13 @@
         ii = ind+m
         mnum=molnum(i)
         iti=itype(i,mnum)
-        massvec(ii)=msc(iabs(iti),mnum)
+        if (mnum.eq.5) then
+        mscab=0.0
+        else
+        mscab=msc(iabs(iti),mnum)
+        endif
+        massvec(ii)=mscab
+
         if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum) .and. mnum.ne.5) then
           ind1=ind1+1
           ii1= ind1+m1