X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fenergy_p_new_barrier.F;h=87555f96f89d161323ad2cb40c7da020f3e1ec01;hb=2958e81241cbcb0ccf96891c89c753689ac6f709;hp=851ad39228d816dd55e058b3235357e0aa363c6f;hpb=4866688f3356059b9a933ec2be6da24e6784fa9f;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 851ad39..87555f9 100644 --- a/source/unres/src_MD/energy_p_new_barrier.F +++ b/source/unres/src_MD/energy_p_new_barrier.F @@ -1065,7 +1065,7 @@ C------------------------------------------------------------------------ & 'ESC= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC local)'/ & 'ETORS= ',1pE16.6,' WEIGHT=',1pD16.6,' (torsional)'/ & 'ETORSD=',1pE16.6,' WEIGHT=',1pD16.6,' (double torsional)'/ - & 'EHBP= ',1pE16.6,' WEIGHT=',1pD16.6, + & 'EHPB= ',1pE16.6,' WEIGHT=',1pD16.6, & ' (SS bridges & dist. cnstr.)'/ & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ @@ -1137,8 +1137,8 @@ C c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=itype(i) - itypi1=itype(i+1) + itypi=iabs(itype(i)) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1151,7 +1151,7 @@ C cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) + itypj=iabs(itype(j)) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -1314,8 +1314,8 @@ C c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=itype(i) - itypi1=itype(i+1) + itypi=iabs(itype(i)) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1324,7 +1324,7 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) + itypj=iabs(itype(j)) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -1431,8 +1431,8 @@ c else c endif ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) - itypi1=itype(i+1) + itypi=iabs(itype(i)) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1567,8 +1567,8 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.eq.0) lprn=.false. ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) - itypi1=itype(i+1) + itypi=iabs(itype(i)) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1585,7 +1585,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, @@ -1726,8 +1726,8 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.eq.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) - itypi1=itype(i+1) + itypi=iabs(itype(i)) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1742,7 +1742,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) @@ -2046,8 +2046,8 @@ C cd print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=itype(i) - itypi1=itype(i+1) + itypi=iabs(itype(i)) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -2058,7 +2058,7 @@ C cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) + itypj=iabs(itype(j)) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -4059,7 +4059,7 @@ cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) - itypj=itype(j) + itypj=iabs(itype(j)) C Uncomment following three lines for SC-p interactions c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi @@ -4153,7 +4153,7 @@ cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) - itypj=itype(j) + itypj=iabs(itype(j)) C Uncomment following three lines for SC-p interactions c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi @@ -4266,49 +4266,91 @@ C iii and jjj point to the residues for which the distance is assigned. iii=ii jjj=jj endif -cd write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj +c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj, +c & dhpb(i),dhpb1(i),forcon(i) C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. - if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then + if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. iabs(itype(jjj + &)).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij cd write (iout,*) "eij",eij + else if (ii.gt.nres .and. jj.gt.nres) then +c Restraints from contact prediction + dd=dist(ii,jj) + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "beta nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else + dd=dist(ii,jj) + rdis=dd-dhpb(i) +C Get the force constant corresponding to this distance. + waga=forcon(i) +C Calculate the contribution to energy. + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "beta reg",dd,waga*rdis*rdis +C +C Evaluate gradient. +C + fac=waga*rdis/dd + endif + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo + do j=1,3 + ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) + ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) + enddo + do k=1,3 + ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) + enddo else C Calculate the distance between the two points and its difference from the C target distance. - dd=dist(ii,jj) - rdis=dd-dhpb(i) + dd=dist(ii,jj) + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "alph nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else + rdis=dd-dhpb(i) C Get the force constant corresponding to this distance. - waga=forcon(i) + waga=forcon(i) C Calculate the contribution to energy. - ehpb=ehpb+waga*rdis*rdis + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "alpha reg",dd,waga*rdis*rdis C C Evaluate gradient. C - fac=waga*rdis/dd + fac=waga*rdis/dd + endif cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, cd & ' waga=',waga,' fac=',fac - do j=1,3 - ggg(j)=fac*(c(j,jj)-c(j,ii)) - enddo + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3) C If this is a SC-SC distance, we need to calculate the contributions to the C Cartesian gradient in the SC vectors (ghpbx). - if (iii.lt.ii) then + if (iii.lt.ii) then do j=1,3 ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) enddo - endif + endif cgrad do j=iii,jjj-1 cgrad do k=1,3 cgrad ghpbc(k,j)=ghpbc(k,j)+ggg(k) cgrad enddo cgrad enddo - do k=1,3 - ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) - ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) - enddo + do k=1,3 + ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) + enddo endif enddo ehpb=0.5D0*ehpb @@ -4333,7 +4375,7 @@ C include 'COMMON.VAR' include 'COMMON.IOUNITS' double precision erij(3),dcosom1(3),dcosom2(3),gg(3) - itypi=itype(i) + itypi=iabs(itype(i)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -4342,7 +4384,7 @@ C dzi=dc_norm(3,nres+i) c dsci_inv=dsc_inv(itypi) dsci_inv=vbld_inv(nres+i) - itypj=itype(j) + itypj=iabs(itype(j)) c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(nres+j) xj=c(1,nres+j)-xi @@ -4436,7 +4478,7 @@ c c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included c do i=ibond_start,ibond_end - iti=itype(i) + iti=iabs(itype(i)) if (iti.ne.10) then nbi=nbondterm(iti) if (nbi.eq.1) then @@ -4511,7 +4553,7 @@ c write (*,'(a,i2)') 'EBEND ICG=',icg do i=ithet_start,ithet_end C Zero the energy function and its derivative at 0 or pi. call splinthet(theta(i),0.5d0*delta,ss,ssd) - it=itype(i-1) + it=iabs(itype(i-1)) if (i.gt.3) then #ifdef OSF phii=phi(i) @@ -4725,7 +4767,7 @@ C dephii=0.0d0 dephii1=0.0d0 theti2=0.5d0*theta(i) - ityp2=ithetyp(itype(i-1)) + ityp2=ithetyp(iabs(itype(i-1))) do k=1,nntheterm coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) @@ -4737,7 +4779,7 @@ C #else phii=phi(i) #endif - ityp1=ithetyp(itype(i-2)) + ityp1=ithetyp(iabs(itype(i-2))) do k=1,nsingle cosph1(k)=dcos(k*phii) sinph1(k)=dsin(k*phii) @@ -4758,7 +4800,7 @@ C #else phii1=phi(i+1) #endif - ityp3=ithetyp(itype(i)) + ityp3=ithetyp(iabs(itype(i))) do k=1,nsingle cosph2(k)=dcos(k*phii1) sinph2(k)=dsin(k*phii1) @@ -4909,7 +4951,7 @@ c write (iout,'(a)') 'ESC' do i=loc_start,loc_end it=itype(i) if (it.eq.10) goto 1 - nlobit=nlob(it) + nlobit=nlob(iabs(it)) c print *,'i=',i,' it=',it,' nlobit=',nlobit c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad theti=theta(i+1)-pipol @@ -5066,11 +5108,11 @@ C Compute the contribution to SC energy and derivatives do j=1,nlobit #ifdef OSF - adexp=bsc(j,it)-0.5D0*contr(j,iii)+emin + adexp=bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin if(adexp.ne.adexp) adexp=1.0 expfac=dexp(adexp) #else - expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin) + expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin) #endif cd print *,'j=',j,' expfac=',expfac escloc_i=escloc_i+expfac @@ -5152,7 +5194,7 @@ C Compute the contribution to SC energy and derivatives dersc12=0.0d0 do j=1,nlobit - expfac=dexp(bsc(j,it)-0.5D0*contr(j)+emin) + expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j)+emin) escloc_i=escloc_i+expfac do k=1,2 dersc(k)=dersc(k)+Ax(k,j)*expfac @@ -5731,12 +5773,17 @@ c lprn=.true. etors_ii=0.0D0 itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) + if (iabs(itype(i)).eq.20) then + iblock=2 + else + iblock=1 + endif phii=phi(i) gloci=0.0D0 C Regular cosine and sine terms - do j=1,nterm(itori,itori1) - v1ij=v1(j,itori,itori1) - v2ij=v2(j,itori,itori1) + do j=1,nterm(itori,itori1,iblock) + v1ij=v1(j,itori,itori1,iblock) + v2ij=v2(j,itori,itori1,iblock) cosphi=dcos(j*phii) sinphi=dsin(j*phii) etors=etors+v1ij*cosphi+v2ij*sinphi @@ -5751,7 +5798,7 @@ C [v2 cos(phi/2)+v3 sin(phi/2)]^2 + 1 C cosphi=dcos(0.5d0*phii) sinphi=dsin(0.5d0*phii) - do j=1,nlor(itori,itori1) + do j=1,nlor(itori,itori1,iblock) vl1ij=vlor1(j,itori,itori1) vl2ij=vlor2(j,itori,itori1) vl3ij=vlor3(j,itori,itori1) @@ -5764,13 +5811,14 @@ C gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom enddo C Subtract the constant term - etors=etors-v0(itori,itori1) + etors=etors-v0(itori,itori1,iblock) if (energy_dec) write (iout,'(a6,i5,0pf7.3)') - & 'etor',i,etors_ii-v0(itori,itori1) + & 'etor',i,etors_ii-v0(itori,itori1,iblock) if (lprn) & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1, - & (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6) + & (v1(j,itori,itori1,iblock),j=1,6), + & (v2(j,itori,itori1,iblock),j=1,6) gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) enddo @@ -8053,9 +8101,9 @@ C C Parallel Antiparallel C C o o -C /l\ /j\ -C / \ / \ -C /| o | | o |\ +C /l\ /j\ +C / \ / \ +C /| o | | o |\ C \ j|/k\| / \ |/k\|l / C \ / \ / \ / \ / C o o o o @@ -8154,18 +8202,18 @@ c---------------------------------------------------------------------------- logical lprn common /kutas/ lprn CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C -C Parallel Antiparallel -C -C o o -C \ /l\ /j\ / -C \ / \ / \ / -C o| o | | o |o -C \ j|/k\| \ |/k\|l -C \ / \ \ / \ -C o o -C i i -C +C C +C Parallel Antiparallel C +C C +C o o C +C \ /l\ /j\ / C +C \ / \ / \ / C +C o| o | | o |o C +C \ j|/k\| \ |/k\|l C +C \ / \ \ / \ C +C o o C +C i i C +C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cd write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l C AL 7/4/01 s1 would occur in the sixth-order moment, @@ -8336,18 +8384,18 @@ c---------------------------------------------------------------------------- double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2) logical swap CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C -C Parallel Antiparallel -C -C o o -C /l\ / \ /j\ -C / \ / \ / \ -C /| o |o o| o |\ -C j|/k\| / |/k\|l / -C / \ / / \ / -C / o / o -C i i -C +C C +C Parallel Antiparallel C +C C +C o o C +C /l\ / \ /j\ C +C / \ / \ / \ C +C /| o |o o| o |\ C +C j|/k\| / |/k\|l / C +C / \ / / \ / C +C / o / o C +C i i C +C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 4/7/01 AL Component s1 was removed, because it pertains to the respective @@ -8453,18 +8501,18 @@ c---------------------------------------------------------------------------- & auxvec1(2),auxmat1(2,2) logical swap CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C -C Parallel Antiparallel -C -C o o -C /l\ / \ /j\ -C / \ / \ / \ -C /| o |o o| o |\ -C \ j|/k\| \ |/k\|l -C \ / \ \ / \ -C o \ o \ -C i i -C +C C +C Parallel Antiparallel C +C C +C o o C +C /l\ / \ /j\ C +C / \ / \ / \ C +C /| o |o o| o |\ C +C \ j|/k\| \ |/k\|l C +C \ / \ \ / \ C +C o \ o \ C +C i i C +C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 4/7/01 AL Component s1 was removed, because it pertains to the respective