From 58b37a8c9dff4d0d448f99de2ea09225f105faa0 Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Tue, 22 Mar 2016 04:48:24 +0100 Subject: [PATCH] energy_dec for esccor and some other debug printout --- source/unres/src_MD-M/energy_p_new_barrier.F | 7 ++++++- source/unres/src_MD-M/geomout.F | 2 +- source/unres/src_MD-M/intcor.f | 7 +++++-- source/unres/src_MD-M/lagrangian_lesyng.F | 1 + source/unres/src_MD-M/unres.F | 1 + source/unres/src_MD/energy_p_new_barrier.F | 6 +++++- source/unres/src_MD/lagrangian_lesyng.F | 1 + source/unres/src_MD/unres.F | 1 + 8 files changed, 21 insertions(+), 5 deletions(-) diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 9ea982a..2bea892 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -5747,12 +5747,13 @@ c write (iout,*) "EBACK_SC_COR",itau_start,itau_end esccor=0.0D0 do i=itau_start,itau_end if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle - esccor_ii=0.0D0 + 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 + esccor_ii=0.0D0 cc Added 09 May 2012 (Adasko) cc Intertyp means interaction type of backbone mainchain correlation: c 1 = SC...Ca...Ca...Ca @@ -5776,9 +5777,13 @@ c 3 = SC...Ca...Ca...SCi v2ij=v2sccor(j,intertyp,isccori,isccori1) cosphi=dcos(j*tauangle(intertyp,i)) sinphi=dsin(j*tauangle(intertyp,i)) + if (energy_dec) esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi esccor=esccor+v1ij*cosphi+v2ij*sinphi gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) enddo + if (energy_dec) write (iout,'(a6,i5,i2,0pf7.3)') + & 'esccor',i,intertyp,esccor_ii +cd write (iout,*) "tau ",i,intertyp,tauangle(intertyp,i)*RAD2DEG c write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci if (lprn) diff --git a/source/unres/src_MD-M/geomout.F b/source/unres/src_MD-M/geomout.F index d21b400..de09e9f 100644 --- a/source/unres/src_MD-M/geomout.F +++ b/source/unres/src_MD-M/geomout.F @@ -9,7 +9,7 @@ include 'COMMON.SBRIDGE' include 'COMMON.DISTFIT' include 'COMMON.MD' - character*50 tytul + character*(*) tytul character*1 chainid(10) /'A','B','C','D','E','F','G','H','I','J'/ dimension ica(maxres) write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot diff --git a/source/unres/src_MD-M/intcor.f b/source/unres/src_MD-M/intcor.f index a3cd5d0..04edbfd 100644 --- a/source/unres/src_MD-M/intcor.f +++ b/source/unres/src_MD-M/intcor.f @@ -32,6 +32,7 @@ c include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' x12=c(1,i1)-c(1,i2) x23=c(1,i3)-c(1,i2) x34=c(1,i4)-c(1,i3) @@ -57,8 +58,6 @@ cd print '(2i3,3f10.5)',i3,i4,x34,y34,z34 if (dabs(scalar).gt.1.0D0) &scalar=0.99999999999999D0*scalar/dabs(scalar) angle=dacos(scalar) -cd print '(2i4,10f7.3)',i2,i3,vx,vy,vz,wx,wy,wz,vnorm,wnorm, -cd &scalar,angle else angle=pi endif @@ -69,6 +68,10 @@ c if (angle.le.0.0D0) angle=pi+angle scalar=tx*x23+ty*y23+tz*z23 if (scalar.lt.0.0D0) angle=-angle beta=angle +cd if ((vnorm.lt.0.01.or.wnorm.lt.0.01).and.angle.ne.pi +cd & .and.i3.ne.maxres*2) +cd & write(iout,'(a12,4i4,6f7.3,4f10.5)') 'beta warning',i1,i2,i3,i4, +cd & vx,vy,vz,wx,wy,wz,vnorm,wnorm,scalar,angle*RAD2DEG return end C diff --git a/source/unres/src_MD-M/lagrangian_lesyng.F b/source/unres/src_MD-M/lagrangian_lesyng.F index 3398091..de3858a 100644 --- a/source/unres/src_MD-M/lagrangian_lesyng.F +++ b/source/unres/src_MD-M/lagrangian_lesyng.F @@ -182,6 +182,7 @@ c sites (dimen1) endif #endif c write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3 + write (iout,*) "The number of degrees of freedom ",dimen3 c Zeroing out A and fricmat do i=1,dimen do j=1,dimen diff --git a/source/unres/src_MD-M/unres.F b/source/unres/src_MD-M/unres.F index 83d9588..2df2658 100644 --- a/source/unres/src_MD-M/unres.F +++ b/source/unres/src_MD-M/unres.F @@ -258,6 +258,7 @@ crc overlap test time1=tcpu() #endif call minim_dc(etot,iretcode,nfun) + if(iretcode.eq.8) call check_ecartint else if (indpdb.ne.0) then call bond_regular diff --git a/source/unres/src_MD/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F index 3a7a1f6..19e3581 100644 --- a/source/unres/src_MD/energy_p_new_barrier.F +++ b/source/unres/src_MD/energy_p_new_barrier.F @@ -5989,7 +5989,7 @@ c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor esccor=0.0D0 do i=itau_start,itau_end C do i=42,42 - esccor_ii=0.0D0 + if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle isccori=isccortyp(itype(i-2)) isccori1=isccortyp(itype(i-1)) @@ -6003,6 +6003,7 @@ c(see comment below) C print *,i,tauangle(1,i) do intertyp=1,3 !intertyp + esccor_ii=0.0D0 cc Added 09 May 2012 (Adasko) cc Intertyp means interaction type of backbone mainchain correlation: c 1 = SC...Ca...Ca...Ca @@ -6024,9 +6025,12 @@ c 3 = SC...Ca...Ca...SCi v2ij=v2sccor(j,intertyp,isccori,isccori1) cosphi=dcos(j*tauangle(intertyp,i)) sinphi=dsin(j*tauangle(intertyp,i)) + if (energy_dec) esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi esccor=esccor+v1ij*cosphi+v2ij*sinphi gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) enddo + if (energy_dec) write (iout,'(a6,i5,i2,0pf7.3)') + & 'esccor',i,intertyp,esccor_ii C print *,i,tauangle(1,i),gloci gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci c write (iout,*) "WTF",intertyp,i,itype(i),v1ij*cosphi+v2ij*sinphi diff --git a/source/unres/src_MD/lagrangian_lesyng.F b/source/unres/src_MD/lagrangian_lesyng.F index 775142b..0c188a9 100644 --- a/source/unres/src_MD/lagrangian_lesyng.F +++ b/source/unres/src_MD/lagrangian_lesyng.F @@ -184,6 +184,7 @@ c sites (dimen1) endif #endif c write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3 + write (iout,*) "The number of degrees of freedom ",dimen3 c Zeroing out A and fricmat do i=1,dimen do j=1,dimen diff --git a/source/unres/src_MD/unres.F b/source/unres/src_MD/unres.F index a63047f..722a3aa 100644 --- a/source/unres/src_MD/unres.F +++ b/source/unres/src_MD/unres.F @@ -248,6 +248,7 @@ crc overlap test time1=tcpu() #endif call minim_dc(etot,iretcode,nfun) + if(iretcode.eq.8) call check_ecartint else if (indpdb.ne.0) then call bond_regular -- 1.7.9.5