X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Fenergy_p_new_barrier.F;h=5e90c170ab681644e245cd3200134b4c4395d77e;hb=34d3ad3987785642be58fb2f26557d3314215577;hp=c9140a6ccafcd091dd375babb7a9bdedcb02bd6a;hpb=f690e8b70bab14132839afebf080d4a28363b226;p=unres.git diff --git a/source/unres/src_MD-M-newcorr/energy_p_new_barrier.F b/source/unres/src_MD-M-newcorr/energy_p_new_barrier.F index c9140a6..5e90c17 100644 --- a/source/unres/src_MD-M-newcorr/energy_p_new_barrier.F +++ b/source/unres/src_MD-M-newcorr/energy_p_new_barrier.F @@ -121,6 +121,10 @@ C C Calculate electrostatic (H-bonding) energy of the main chain. C 107 continue +cmc +cmc Sep-06: egb takes care of dynamic ss bonds too +cmc +C if (dyn_ss) call dyn_set_nss c print *,"Processor",myrank," computed USCSC" #ifdef TIMING time01=MPI_Wtime() @@ -300,6 +304,7 @@ c Here are the energies showed per procesor if the are more processors c per molecule then we sum it up in sum_energy subroutine c print *," Processor",myrank," calls SUM_ENERGY" call sum_energy(energia,.true.) + if (dyn_ss) call dyn_set_nss c print *," Processor",myrank," left SUM_ENERGY" #ifdef TIMING time_sumene=time_sumene+MPI_Wtime()-time00 @@ -710,6 +715,7 @@ c enddo do i=1,4*nres glocbuf(i)=gloc(i,icg) enddo +#define DEBUG #ifdef DEBUG write (iout,*) "gloc_sc before reduce" do i=1,nres @@ -718,6 +724,7 @@ c enddo enddo enddo #endif +#undef DEBUG do i=1,nres do j=1,3 gloc_scbuf(j,i)=gloc_sc(j,i,icg) @@ -737,6 +744,7 @@ c enddo call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres, & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) time_reduce=time_reduce+MPI_Wtime()-time00 +#define DEBUG #ifdef DEBUG write (iout,*) "gloc_sc after reduce" do i=1,nres @@ -745,6 +753,7 @@ c enddo enddo enddo #endif +#undef DEBUG #ifdef DEBUG write (iout,*) "gloc after reduce" do i=1,4*nres @@ -1408,6 +1417,7 @@ C include 'COMMON.IOUNITS' include 'COMMON.CALC' include 'COMMON.CONTROL' + include 'COMMON.SBRIDGE' logical lprn evdw=0.0D0 ccccc energy_dec=.false. @@ -1435,6 +1445,12 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) + IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN + call dyn_ssbond_ene(i,j,evdwij) + evdw=evdw+evdwij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') + & 'evdw',i,j,evdwij,' ss' + ELSE ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -1529,6 +1545,7 @@ C Calculate the radial part of the gradient gg(3)=zj*fac C Calculate angular part of the gradient. call sc_grad + endif ! dyn_ss enddo ! j enddo ! iint enddo ! i @@ -2267,27 +2284,27 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then iti1=ntortyp+1 endif c write(iout,*),i - b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0d0) - & +bnew1(2,1,iti)*dsin(theta(i-1)) - & +bnew1(3,1,iti)*dcos(theta(i-1)/2.0d0) - gtb1(1,i-2)=bnew1(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0 - & +bnew1(2,1,iti)*dcos(theta(i-1)) - & -bnew1(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0 -c & +bnew1(3,1,iti)*dsin(alpha(i))*cos(beta(i)) + b1(1,i-2)=bnew1(1,1,iti)*sin(theta(i-1)/2.0) + & +bnew1(2,1,iti)*sin(theta(i-1)) + & +bnew1(3,1,iti)*cos(theta(i-1)/2.0) + gtb1(1,i-2)=bnew1(1,1,iti)*cos(theta(i-1)/2.0)/2.0 + & +bnew1(2,1,iti)*cos(theta(i-1)) + & -bnew1(3,1,iti)*sin(theta(i-1)/2.0)/2.0 +c & +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i)) c &*(cos(theta(i)/2.0) - b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0d0) - & +bnew2(2,1,iti)*dsin(theta(i-1)) - & +bnew2(3,1,iti)*dcos(theta(i-1)/2.0d0) -c & +bnew2(3,1,iti)*dsin(alpha(i))*dcos(beta(i)) + b2(1,i-2)=bnew2(1,1,iti)*sin(theta(i-1)/2.0) + & +bnew2(2,1,iti)*sin(theta(i-1)) + & +bnew2(3,1,iti)*cos(theta(i-1)/2.0) +c & +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i)) c &*(cos(theta(i)/2.0) - gtb2(1,i-2)=bnew2(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0 - & +bnew2(2,1,iti)*dcos(theta(i-1)) - & -bnew2(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0 + gtb2(1,i-2)=bnew2(1,1,iti)*cos(theta(i-1)/2.0)/2.0 + & +bnew2(2,1,iti)*cos(theta(i-1)) + & -bnew2(3,1,iti)*sin(theta(i-1)/2.0)/2.0 c if (ggb1(1,i).eq.0.0d0) then c write(iout,*) 'i=',i,ggb1(1,i), -c &bnew1(1,1,iti)*dcos(theta(i)/2.0d0)/2.0d0, -c &bnew1(2,1,iti)*dcos(theta(i)), -c &bnew1(3,1,iti)*dsin(theta(i)/2.0d0)/2.0d0 +c &bnew1(1,1,iti)*cos(theta(i)/2.0)/2.0, +c &bnew1(2,1,iti)*cos(theta(i)), +c &bnew1(3,1,iti)*sin(theta(i)/2.0)/2.0 c endif b1(2,i-2)=bnew1(1,2,iti) gtb1(2,i-2)=0.0 @@ -2304,8 +2321,8 @@ c endif c EE(2,2,iti)=0.0d0 c EE(1,2,iti)=0.5d0*eenew(1,iti) c EE(2,1,iti)=0.5d0*eenew(1,iti) -c b1(2,iti)=bnew1(1,2,iti)*dsin(alpha(i))*dsin(beta(i)) -c b2(2,iti)=bnew2(1,2,iti)*dsin(alpha(i))*dsin(beta(i)) +c b1(2,iti)=bnew1(1,2,iti)*sin(alpha(i))*sin(beta(i)) +c b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i)) b1tilde(1,i-2)=b1(1,i-2) b1tilde(2,i-2)=-b1(2,i-2) b2tilde(1,i-2)=b2(1,i-2) @@ -2319,17 +2336,6 @@ c write (iout,*) 'theta=', theta(i-1) do i=3,nres+1 #endif #endif - if (i.gt. nnt+2 .and. i.lt.nct+2) then - iti = itortyp(itype(i-2)) - else - iti=ntortyp+1 - endif -c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then - if (i.gt. nnt+1 .and. i.lt.nct+1) then - iti1 = itortyp(itype(i-1)) - else - iti1=ntortyp+1 - endif if (i .lt. nres+1) then sin1=dsin(phi(i)) cos1=dcos(phi(i)) @@ -2466,12 +2472,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 #ifdef MUOUT - write (iout,'(2hmu,i3,3f8.1,12f10.5)') i-2,rad2deg*theta(i-1), + write (iout,'(2hmu,i3,3f8.1,7f10.5)') i-2,rad2deg*theta(i-1), & rad2deg*theta(i),rad2deg*phi(i),mu(1,i-2),mu(2,i-2), & -b2(1,i-2),b2(2,i-2),b1(1,i-2),b1(2,i-2), & dsqrt(b2(1,i-1)**2+b2(2,i-1)**2) - & +dsqrt(b1(1,i-1)**2+b1(2,i-1)**2), - & ((ee(l,k,i-2),l=1,2),k=1,2),eenew(1,itortyp(iti)) + & +dsqrt(b1(1,i-1)**2+b1(2,i-1)**2) #endif cd write (iout,*) 'mu ',mu(:,i-2) cd write (iout,*) 'mu1',mu1(:,i-2) @@ -4288,10 +4293,20 @@ C iii and jjj point to the residues for which the distance is assigned. cd write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj 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. iabs(itype(iii)).eq.1 .and. - & iabs(itype(jjj)).eq.1) then - call ssbond_ene(iii,jjj,eij) +cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then +C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds +c if (.not.dyn_ss .and. i.le.nss) then +C 15/02/13 CC dynamic SSbond +C if (.not.dyn_ss.and. +C & ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then + + if (.not.dyn_ss .and. i.le.nss) then +C 15/02/13 CC dynamic SSbond - additional check + if (ii.gt.nres + & .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then + call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij + endif cd write (iout,*) "eij",eij else C Calculate the distance between the two points and its difference from the @@ -5804,7 +5819,7 @@ c---------------------------------------------------------------------------- logical lprn C Set lprn=.true. for debugging lprn=.false. -c lprn=.true. +c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 @@ -5852,12 +5867,12 @@ C C Subtract the constant term etors=etors-v0(itori,itori1,iblock) if (energy_dec) write (iout,'(a6,i5,0pf7.3)') - & 'etor',i,etors_ii + & 'etor',i,etors_ii-v0(itori,itori1,iblock) if (lprn) - & write (iout,'(2(a3,2x,i3,2x),2i3,f10.2,6f8.3/36x,6f8.3/)') + & 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, - & rad2deg*phii, - & (v1(j,itori,itori1,1),j=1,6),(v2(j,itori,itori1,1),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