X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fenergy_p_new_barrier.F;h=fbc47917f5a150d38f28e714a619a5e2316fe540;hb=ecf5b7759f8fd1e9656b3bbbe8b4f05a81410aa2;hp=0bcd6f3366ef134c468492a68eca8435a1040045;hpb=d5132ba348c64808e25b811a93271b9f01f8edb8;p=unres.git diff --git a/source/unres/src_MD/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F index 0bcd6f3..fbc4791 100644 --- a/source/unres/src_MD/energy_p_new_barrier.F +++ b/source/unres/src_MD/energy_p_new_barrier.F @@ -1679,9 +1679,10 @@ c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 & evdwij endif - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'evdw',i,j,evdwij - + if (energy_dec) then + write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij + call flush(iout) + endif C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 fac=-expon*(e1+evdwij)*rij_shift @@ -3023,6 +3024,9 @@ C C Loop over i,i+2 and i,i+3 pairs of the peptide groups C do i=iturn3_start,iturn3_end +C if (itype(i).eq.21 .or. itype(i+1).eq.21 +C & .or. itype(i+2).eq.21 .or. itype(i+3).eq.21.or.itype(i+4).eq.21) +C & cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -3038,6 +3042,10 @@ C num_cont_hb(i)=num_conti enddo do i=iturn4_start,iturn4_end +C if (itype(i).eq.21 .or. itype(i+1).eq.21 +C & .or. itype(i+2).eq.21 .or. itype(i+3).eq.21.or.itype(i+4).eq.21 +C & .or. itype(i+5).eq.21) +C & cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -3056,6 +3064,8 @@ 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 +C if (itype(i).eq.21 .or. itype(i+1).eq.21 +C &.or.itype(i+2)) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -3068,6 +3078,8 @@ c c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) do j=ielstart(i),ielend(i) +C if (itype(j).eq.21 .or. itype(j+1).eq.21 +C &.or.itype(j+2)) cycle call eelecij(i,j,ees,evdw1,eel_loc) enddo ! j num_cont_hb(i)=num_conti @@ -4482,6 +4494,8 @@ c do i=ibondp_start,ibondp_end diff = vbld(i)-vbldp0 c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff + if (energy_dec) write (iout,'(a7,i5,4f7.3)') + & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff estr=estr+diff*diff do j=1,3 gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i) @@ -4500,6 +4514,12 @@ c diff=vbld(i+nres)-vbldsc0(1,iti) c write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff, c & AKSC(1,iti),AKSC(1,iti)*diff*diff + if (energy_dec) then + write (iout,*) + & "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff, + & AKSC(1,iti),AKSC(1,iti)*diff*diff + call flush(iout) + endif estr=estr+0.5d0*AKSC(1,iti)*diff*diff do j=1,3 gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres) @@ -4790,7 +4810,7 @@ C sinkt(k)=dsin(k*theti2) enddo C if (i.gt.3) then - if (i.gt.3 .and. itype(i-3).ne.ntyp1) then + if (i.gt.3 .and. itype(imax0(i-3,1)).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -4931,12 +4951,12 @@ C if (i.gt.3) then enddo enddo 10 continue -c lprn1=.true. if (lprn1) write (iout,'(a4,i2,3f8.1,9h ethetai ,f10.5)') & 'ebe', i,theta(i)*rad2deg,phii*rad2deg, & phii1*rad2deg,ethetai -c lprn1=.false. etheta=etheta+ethetai + if (energy_dec) write (iout,'(a6,i5,0pf7.3)') + & 'ebend',i,ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 gloc(nphi+i-2,icg)=wang*dethetai @@ -5878,12 +5898,14 @@ C 6/23/01 Compute double torsional energy include 'COMMON.IOUNITS' include 'COMMON.FFIELD' include 'COMMON.TORCNSTR' + include 'COMMON.CONTROL' logical lprn C Set lprn=.true. for debugging lprn=.false. c lprn=.true. etors_d=0.0D0 do i=iphid_start,iphid_end + etors_d_ii=0.0D0 itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) @@ -5902,6 +5924,8 @@ c lprn=.true. sinphi2=dsin(j*phii1) etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+ & v2cij*cosphi2+v2sij*sinphi2 + if (energy_dec) etors_d_ii=etors_d_ii+ + & v1cij*cosphi1+v1sij*sinphi1+v2cij*cosphi2+v2sij*sinphi2 gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1) gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2) enddo @@ -5917,12 +5941,17 @@ c lprn=.true. sinphi1m2=dsin(l*phii-(k-l)*phii1) etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+ & v1sdij*sinphi1p2+v2sdij*sinphi1m2 + if (energy_dec) etors_d_ii=etors_d_ii+ + & v1cdij*cosphi1p2+v2cdij*cosphi1m2+ + & v1sdij*sinphi1p2+v2sdij*sinphi1m2 gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2 & -v1cdij*sinphi1p2-v2cdij*sinphi1m2) gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2 & -v1cdij*sinphi1p2+v2cdij*sinphi1m2) enddo enddo + if (energy_dec) write (iout,'(a6,i5,0pf7.3)') + & 'etor_d',i,etors_d_ii gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1 gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2 c write (iout,*) "gloci", gloc(i-3,icg) @@ -5959,19 +5988,22 @@ c lprn=.true. 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)) phii=phi(i) + cccc Added 9 May 2012 cc Tauangle is torsional engle depending on the value of first digit c(see comment below) cc Omicron is flat angle depending on the value of first digit c(see comment below) - +C print *,i,tauangle(1,i) - do intertyp=1,3 !intertyp +c do intertyp=1,3 !intertyp + do intertyp=2,2 !intertyp cc Added 09 May 2012 (Adasko) cc Intertyp means interaction type of backbone mainchain correlation: c 1 = SC...Ca...Ca...Ca @@ -5996,6 +6028,7 @@ c 3 = SC...Ca...Ca...SCi esccor=esccor+v1ij*cosphi+v2ij*sinphi gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) enddo +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 c &gloc_sc(intertyp,i-3,icg) @@ -6008,7 +6041,8 @@ c &gloc_sc(intertyp,i-3,icg) enddo !intertyp enddo c do i=1,nres -c write (iout,*) "W@T@F", gloc_sc(1,i,icg),gloc(i,icg) +c write (iout,*) "W@T@F", gloc_sc(1,i,icg),gloc_sc(2,i,icg), +c & gloc_sc(3,i,icg) c enddo return end