X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=6fb9591c6bfd3fa8943c79772786e4ce69f9f88d;hb=e93f349536ec939eb4399c65cec887af423a4eec;hp=79cd1b9815dbf5a43525e4da5f3d5122684b9c2c;hpb=6304fe426dd266ba99503f2cc4a5b9d396511b2f;p=unres.git 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 79cd1b9..6fb9591 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -2860,7 +2860,8 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then cd write (iout,*) '*******i',i,' iti1',iti cd write (iout,*) 'b1',b1(:,iti) cd write (iout,*) 'b2',b2(:,iti) -cd write (iout,*) 'Ug',Ug(:,:,i-2) +cd write (iout,*) "phi(",i,")=",phi(i)," sin1",sin1," cos1",cos1 +cd write (iout,*) 'Ug',Ug(:,:,i-2) c if (i .gt. iatel_s+2) then if (i .gt. nnt+2) then call matvec2(Ug(1,1,i-2),b2(1,i-2),Ub2(1,i-2)) @@ -2915,7 +2916,11 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1) enddo C write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2) -c write (iout,*) 'mu ',mu(:,i-2),i-2 +cd write (iout,*) 'mu ',mu(:,i-2),i-2 +cd write (iout,*) 'b1 ',b1(:,i-1),i-2 +cd write (iout,*) 'Ub2 ',Ub2(:,i-2),i-2 +cd write (iout,*) 'Ug ',Ug(:,:,i-2),i-2 +cd write (iout,*) 'b2 ',b2(:,i-2),i-2 cd write (iout,*) 'mu1',mu1(:,i-2) cd write (iout,*) 'mu2',mu2(:,i-2) if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) @@ -3324,21 +3329,21 @@ C Loop over i,i+2 and i,i+3 pairs of the peptide groups C C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition do i=iturn3_start,iturn3_end - if (i.le.1) cycle +CAna if (i.le.1) cycle C write(iout,*) "tu jest i",i if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((i+4).gt.nres) - & .or.((i-1).le.0) +CAna & .or.((i+4).gt.nres) +CAna & .or.((i-1).le.0) C end of changes by Ana & .or. itype(i+2).eq.ntyp1 & .or. itype(i+3).eq.ntyp1) cycle - if(i.gt.1)then - if(itype(i-1).eq.ntyp1)cycle - end if - if(i.LT.nres-3)then - if (itype(i+4).eq.ntyp1) cycle - end if +CAna if(i.gt.1)then +CAna if(itype(i-1).eq.ntyp1)cycle +CAna end if +CAna if(i.LT.nres-3)then +CAna if (itype(i+4).eq.ntyp1) cycle +CAna end if dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -3360,17 +3365,17 @@ C end of changes by Ana num_cont_hb(i)=num_conti enddo do i=iturn4_start,iturn4_end - if (i.le.1) cycle +cAna if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((i+5).gt.nres) - & .or.((i-1).le.0) +cAna & .or.((i+5).gt.nres) +cAna & .or.((i-1).le.0) C end of changes suggested by Ana & .or. itype(i+3).eq.ntyp1 & .or. itype(i+4).eq.ntyp1 - & .or. itype(i+5).eq.ntyp1 - & .or. itype(i).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +cAna & .or. itype(i+5).eq.ntyp1 +cAna & .or. itype(i).eq.ntyp1 +cAna & .or. itype(i-1).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -3428,14 +3433,14 @@ c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c do i=iatel_s,iatel_e - if (i.le.1) cycle +cAna if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((i+2).gt.nres) - & .or.((i-1).le.0) +cAna & .or.((i+2).gt.nres) +cAna & .or.((i-1).le.0) C end of changes by Ana - & .or. itype(i+2).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +cAna & .or. itype(i+2).eq.ntyp1 +cAna & .or. itype(i-1).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -3486,14 +3491,14 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) do j=ielstart(i),ielend(i) C write (iout,*) i,j - if (j.le.1) cycle +cAna if (j.le.1) cycle if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((j+2).gt.nres) - & .or.((j-1).le.0) +cAna & .or.((j+2).gt.nres) +cAna & .or.((j-1).le.0) C end of changes by Ana - & .or.itype(j+2).eq.ntyp1 - & .or.itype(j-1).eq.ntyp1 +cAna & .or.itype(j+2).eq.ntyp1 +cAna & .or.itype(j-1).eq.ntyp1 &) cycle call eelecij(i,j,ees,evdw1,eel_loc) enddo ! j @@ -5684,7 +5689,7 @@ C logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 do i=ithet_start,ithet_end - if (i.eq.2) cycle +c if (i.eq.2) cycle c print *,i,itype(i-1),itype(i),itype(i-2) if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1) & .or.(itype(i).eq.ntyp1)) cycle @@ -5701,7 +5706,7 @@ C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(i-3).ne.ntyp1) then + if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -5716,7 +5721,7 @@ C propagation of chirality for glycine type enddo else phii=0.0d0 - ityp1=nthetyp+1 + ityp1=ithetyp(itype(i-2)) do k=1,nsingle cosph1(k)=0.0d0 sinph1(k)=0.0d0 @@ -5737,7 +5742,7 @@ C propagation of chirality for glycine type enddo else phii1=0.0d0 - ityp3=nthetyp+1 + ityp3=ithetyp(itype(i)) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@ -6295,6 +6300,8 @@ c & sumene4, c & dscp1,dscp2,sumene c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1)) escloc = escloc + sumene + if (energy_dec) write (iout,'(a6,i5,0pf7.3)') + & 'escloc',i,sumene c write (2,*) "i",i," escloc",sumene,escloc,it,itype(i) c & ,zz,xx,yy c#define DEBUG @@ -6896,7 +6903,18 @@ c endif enddo - min_odl=minval(distancek) +c min_odl=minval(distancek) + do kk=1,constr_homology + if(l_homo(kk,ii)) then + min_odl=distancek(kk) + exit + endif + enddo + do kk=1,constr_homology + if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) + & min_odl=distancek(kk) + enddo + c write (iout,* )"min_odl",min_odl #ifdef DEBUG write (iout,*) "ij dij",i,j,dij @@ -7049,11 +7067,12 @@ c write (iout,*) idihconstr_start_homo,idihconstr_end_homo do i=idihconstr_start_homo,idihconstr_end_homo kat2=0.0d0 c betai=beta(i,i+1,i+2,i+3) - betai = phi(i+3) + betai = phi(i) c write (iout,*) "betai =",betai do k=1,constr_homology dih_diff(k)=pinorm(dih(k,i)-betai) -c write (iout,*) "dih_diff(",k,") =",dih_diff(k) +cd write (iout,'(a8,2i4,2f15.8)') "dih_diff",i,k,dih_diff(k) +cd & ,sigma_dih(k,i) c if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)= c & -(6.28318-dih_diff(i,k)) c if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)= @@ -7096,7 +7115,7 @@ c grad_dih3=sum_sgdih/sum_gdih c write(iout,*)i,k,gdih,sgdih,beta(i+1,i+2,i+3,i+4),grad_dih3 ccc write(iout,747) "GRAD_KAT_1", i, nphi, icg, grad_dih3, ccc & gloc(nphi+i-3,icg) - gloc(i,icg)=gloc(i,icg)+grad_dih3 + gloc(i-3,icg)=gloc(i-3,icg)+grad_dih3 c if (i.eq.25) then c write(iout,*) "i",i,"icg",icg,"gloc(",i,icg,")",gloc(i,icg) c endif @@ -7161,6 +7180,9 @@ c c dtheta_i=theta(j)-thetaref(j,iref) c dtheta_i=thetaref(k,i)-theta(i) ! original form without indexing theta_diff(k)=thetatpl(k,i)-theta(i) +cd write (iout,'(a8,2i4,2f15.8)') "theta_diff",i,k,theta_diff(k) +cd & ,sigma_theta(k,i) + c utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument c utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta? @@ -7231,7 +7253,7 @@ c Original sign inverted for calc of gradient dyy=-yytpl(k,i)+yytab(i) ! ibid y dzz=-zztpl(k,i)+zztab(i) ! ibid z c write(iout,*) "dxx, dyy, dzz" -c write(iout,'(2i5,3f8.2)') k,i,dxx,dyy,dzz +cd write(iout,'(2i5,4f8.2)') k,i,dxx,dyy,dzz,sigma_d(k,i) c usc_diff_i=-0.5d0*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d rmvd from Gaussian argument c usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d? @@ -7535,12 +7557,12 @@ 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 @@ -7564,9 +7586,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 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)