X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;fp=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=747b13261294fcdaf24d1f64271bc3ff048ded44;hb=e4035e50115ab9c0e65b0445aed96096a5cb86d5;hp=00862f24900104f5e82d9398b6736fa125d10a1d;hpb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;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 00862f2..747b132 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -130,6 +130,8 @@ cmc c if (dyn_ss) call dyn_set_nss c print *,"Processor",myrank," computed USCSC" +c write (iout,*) "SCSC computed OK" +c call flush_(iout) #ifdef TIMING time01=MPI_Wtime() #endif @@ -171,6 +173,8 @@ c print *,"Processor",myrank," left VEC_AND_DERIV" c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, c & eello_turn4) endif +c write (iout,*) "eelec computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed UELEC" C C Calculate excluded-volume interaction energy between peptide groups @@ -187,10 +191,14 @@ C c write (iout,*) "Soft-sphere SCP potential" call escp_soft_sphere(evdw2,evdw2_14) endif +c write (iout,*) "escp computed OK" +c call flush_(iout) c c Calculate the bond-stretching energy c call ebond(estr) +c write (iout,*) "ebond computed OK" +c call flush_(iout) C C Calculate the disulfide-bridge and other energy and the contributions C from other distance constraints. @@ -206,12 +214,16 @@ C ebe=0 ethetacnstr=0 endif +c write (iout,*) "ebend computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed UB" C C Calculate the SC local energy. C C print *,"TU DOCHODZE?" call esc(escloc) +c write (iout,*) "esc computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed USC" C C Calculate the virtual-bond torsional energy. @@ -223,6 +235,8 @@ cd print *,'nterm=',nterm etors=0 edihcnstr=0 endif +c write (iout,*) "etor computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed Utor" C C 6/23/01 Calculate double-torsional energy @@ -232,6 +246,8 @@ C else etors_d=0 endif +c write (iout,*) "etor_d computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed Utord" C C 21/5/07 Calculate local sicdechain correlation energy @@ -241,6 +257,8 @@ C else esccor=0.0d0 endif +c write (iout,*) "eback_sc_corr computed OK" +c call flush_(iout) C print *,"PRZED MULIt" c print *,"Processor",myrank," computed Usccorr" C @@ -250,9 +268,13 @@ C n_corr1=0 if ((wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 & .or. wturn6.gt.0.0d0) .and. ipot.lt.6) then +c write (iout,*) "Calling multibody_eello" +c call flush_(iout) call multibody_eello(ecorr,ecorr5,ecorr6,eturn6,n_corr,n_corr1) -cd write(2,*)'multibody_eello n_corr=',n_corr,' n_corr1=',n_corr1, -cd &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6 +c write(iout,*) +c & 'multibody_eello n_corr=',n_corr,' n_corr1=',n_corr1, +c & " ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6 +c call flush_(iout) else ecorr=0.0d0 ecorr5=0.0d0 @@ -260,9 +282,14 @@ cd &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6 eturn6=0.0d0 endif if ((wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) .and. ipot.lt.6) then +c write (iout,*) "Calling multibody_gb_ecorr" +c call flush_(iout) call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1) -cd write (iout,*) "multibody_hb ecorr",ecorr +c write (iout,*) "Exited multibody_hb ecorr",ecorr +c call flush_(iout) endif +c write (iout,*) "multibody computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed Ucorr" C C If performing constraint dynamics, call the constraint energy @@ -281,12 +308,15 @@ C print *,"przed lipidami" if (wliptran.gt.0) then call Eliptransfer(eliptran) endif -C print *,"za lipidami" +c write (iout,*) "lipid energy computed OK" +c call flush_(iout) if (AFMlog.gt.0) then call AFMforce(Eafmforce) else if (selfguide.gt.0) then call AFMvel(Eafmforce) endif +c write (iout,*) "AFMforce computed OK" +c call flush_(iout) #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 #endif @@ -335,7 +365,11 @@ 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.) +c write (iout,*) "sum energy OK" +c call flush_(iout) if (dyn_ss) call dyn_set_nss +c write (iout,*) "Exiting energy" +c call flush_(iout) c print *," Processor",myrank," left SUM_ENERGY" #ifdef TIMING time_sumene=time_sumene+MPI_Wtime()-time00 @@ -5638,6 +5672,8 @@ c time12=1.0d0 etheta=0.0D0 c write (*,'(a,i2)') 'EBEND ICG=',icg do i=ithet_start,ithet_end +c write (iout,*) "ebend: i=",i +c call flush_(iout) if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 & .or.itype(i).eq.ntyp1) cycle C Zero the energy function and its derivative at 0 or pi. @@ -5739,9 +5775,12 @@ C Derivatives of the "mean" values in gamma1 and gamma2. if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg) enddo +c write (iout,*) "Exit loop" +c call flush_(iout) ethetacnstr=0.0d0 -C print *,ithetaconstr_start,ithetaconstr_end,"TU" - do i=ithetaconstr_start,ithetaconstr_end +c write (iout,*) ithetaconstr_start,ithetaconstr_end,"TU" +c call flush_(iout) + do i=max0(ithetaconstr_start,1),ithetaconstr_end itheta=itheta_constr(i) thetiii=theta(itheta) difi=pinorm(thetiii-theta_constr0(i)) @@ -5766,6 +5805,8 @@ C print *,ithetaconstr_start,ithetaconstr_end,"TU" & gloc(itheta+nphi-2,icg) endif enddo +c write (iout,*) "Exit ebend" +c call flush_(iout) C Ufff.... We've done all this!!! return @@ -6103,7 +6144,7 @@ c lprn1=.false. C now constrains ethetacnstr=0.0d0 C print *,ithetaconstr_start,ithetaconstr_end,"TU" - do i=ithetaconstr_start,ithetaconstr_end + do i=max0(ithetaconstr_start,1),ithetaconstr_end itheta=itheta_constr(i) thetiii=theta(itheta) difi=pinorm(thetiii-theta_constr0(i)) @@ -6986,7 +7027,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 C ANY TWO ARE DUMMY ATOMS in row CYCLE @@ -7379,7 +7420,7 @@ C Set lprn=.true. for debugging if (lprn) then write (iout,'(a)') 'Contact function values before RECEIVE:' do i=nnt,nct-2 - write (iout,'(2i3,50(1x,i2,f5.2))') + write (iout,'(2i3,50(1x,i3,f5.2))') & i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), & j=1,num_cont_hb(i)) enddo @@ -7618,6 +7659,7 @@ C Calculate the local-electrostatic correlation terms jp1=iabs(j1) c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, c & ' jj=',jj,' kk=',kk +c call flush(iout) if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0 & .or. j.lt.0 .and. j1.gt.0) .and. & (jp1.eq.jp+1 .or. jp1.eq.jp-1)) then @@ -7635,8 +7677,9 @@ c ecorr=ecorr+ehbcorr(i,j,i+1,j,jj,kk,0.60D0,-0.40D0) enddo ! kk do kk=1,num_conti j1=jcont_hb(kk,i) -c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, -c & ' jj=',jj,' kk=',kk +c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, +c & ' jj=',jj,' kk=',kk +c call flush(iout) if (j1.eq.j+1) then C Contacts I-J and (I+1)-J occur simultaneously. C The system loses extra energy. @@ -7730,6 +7773,7 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding include 'COMMON.CONTACTS' include 'COMMON.CHAIN' include 'COMMON.CONTROL' + include 'COMMON.TORSION' double precision gx(3),gx1(3) integer num_cont_hb_old(maxres) logical lprn,ldone @@ -7738,6 +7782,8 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding C Set lprn=.true. for debugging lprn=.false. eturn6=0.0d0 +c write (iout,*) "MULTIBODY_EELLO" +c call flush(iout) #ifdef MPI do i=1,nres num_cont_hb_old(i)=num_cont_hb(i) @@ -7748,12 +7794,12 @@ C Set lprn=.true. for debugging if (lprn) then write (iout,'(a)') 'Contact function values before RECEIVE:' do i=nnt,nct-2 - write (iout,'(2i3,50(1x,i2,f5.2))') + write (iout,'(2i3,50(1x,i3,f5.2))') & i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), & j=1,num_cont_hb(i)) enddo + call flush(iout) endif - call flush(iout) do i=1,ntask_cont_from ncont_recv(i)=0 enddo @@ -7955,10 +8001,15 @@ c call flush(iout) if (lprn) then write (iout,'(a)') 'Contact function values:' do i=nnt,nct-2 - write (iout,'(2i3,50(1x,i2,5f6.3))') + write (iout,'(2i3,50(1x,i3,5f6.3))') & i,num_cont_hb(i),(jcont_hb(j,i),d_cont(j,i), & ((a_chuj(ll,kk,j,i),ll=1,2),kk=1,2),j=1,num_cont_hb(i)) enddo + write (iout,*) "itortyp" + do i=1,nres + write (iout,*) i,itype(i),itortyp(itype(i)) + enddo + call flush(iout) endif ecorr=0.0D0 ecorr5=0.0d0 @@ -7999,8 +8050,11 @@ c write (iout,*) "corr loop i",i do kk=1,num_conti1 j1=jcont_hb(kk,i1) jp1=iabs(j1) -c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, -c & ' jj=',jj,' kk=',kk + if (lprn) then + write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, + & ' jj=',jj,' kk=',kk + call flush(iout) + endif c if (j1.eq.j+1 .or. j1.eq.j-1) then if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0 & .or. j.lt.0 .and. j1.gt.0) .and. @@ -8040,6 +8094,12 @@ c do iii=1,nres c write (iout,'(i5,3f10.5)') c & iii,(gradcorr5(jjj,iii),jjj=1,3) c enddo +c write (iout,*) "ecorr4" +c call flush(iout) +c write (iout,*) "eello5:",i,jp,i+1,jp1,jj,kk, +c & itype(jp),itype(i+1),itype(jp1), +c & itortyp(itype(jp)),itortyp(itype(i+1)),itortyp(itype(jp1)) +c call flush(iout) if (wcorr5.gt.0.0d0) & ecorr5=ecorr5+eello5(i,jp,i+1,jp1,jj,kk) c write (iout,*) "gradcorr5 after eello5" @@ -8052,12 +8112,16 @@ c enddo 2 'ecorr5',i,j,i+1,j1,eello5(i,jp,i+1,jp1,jj,kk) cd write(2,*)'wcorr6',wcorr6,' wturn6',wturn6 cd write(2,*)'ijkl',i,jp,i+1,jp1 +c write (iout,*) "ecorr5" +c call flush(iout) if (wcorr6.gt.0.0d0 .and. (jp.ne.i+4 .or. jp1.ne.i+3 & .or. wturn6.eq.0.0d0))then cd write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1 ecorr6=ecorr6+eello6(i,jp,i+1,jp1,jj,kk) if (energy_dec) write (iout,'(a6,4i5,0pf7.3)') 1 'ecorr6',i,j,i+1,j1,eello6(i,jp,i+1,jp1,jj,kk) +c write (iout,*) "ecorr6" +c call flush(iout) cd write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5, cd & 'ecorr6=',ecorr6 cd write (iout,'(4e15.5)') sred_geom, @@ -8071,10 +8135,13 @@ cd write (iout,*) '******eturn6: i,j,i+1,j1',i,jip,i+1,jp1 if (energy_dec) write (iout,'(a6,4i5,0pf7.3)') 1 'eturn6',i,j,i+1,j1,eello_turn6(i,jj,kk) cd write (2,*) 'multibody_eello:eturn6',eturn6 +c write (iout,*) "ecorr4" +c call flush(iout) endif ENDIF 1111 continue endif + if (energy_dec) call flush(iout) enddo ! kk enddo ! jj enddo ! i @@ -8837,9 +8904,9 @@ cd endif cd write (iout,*) cd & 'EELLO5: Contacts have occurred for peptide groups',i,j, cd & ' and',k,l - itk=itortyp(itype(k)) - itl=itortyp(itype(l)) - itj=itortyp(itype(j)) +c itk=itortyp(itype(k)) +c itl=itortyp(itype(l)) +c itj=itortyp(itype(j)) eello5_1=0.0d0 eello5_2=0.0d0 eello5_3=0.0d0 @@ -8908,7 +8975,7 @@ C Cartesian gradient c goto 1112 c1111 continue C Contribution from graph II - call transpose2(EE(1,1,itk),auxmat(1,1)) + call transpose2(EE(1,1,k),auxmat(1,1)) call matmat2(auxmat(1,1),AEA(1,1,1),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -8989,7 +9056,7 @@ C Cartesian gradient cd goto 1112 C Contribution from graph IV cd1110 continue - call transpose2(EE(1,1,itl),auxmat(1,1)) + call transpose2(EE(1,1,l),auxmat(1,1)) call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -9062,7 +9129,7 @@ C Cartesian gradient cd goto 1112 C Contribution from graph IV 1110 continue - call transpose2(EE(1,1,itj),auxmat(1,1)) + call transpose2(EE(1,1,j),auxmat(1,1)) call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -9359,7 +9426,7 @@ C o o o o C C i i C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC - itk=itortyp(itype(k)) +c itk=itortyp(itype(k)) s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i)) s2=-scalar2(AEAb2(1,1,imat),Ug2Db1t(1,k)) s3= scalar2(AEAb2(1,1,imat),CUgb2(1,k)) @@ -9649,19 +9716,19 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 4/7/01 AL Component s1 was removed, because it pertains to the respective C energy moment and not to the cluster cumulant. - iti=itortyp(itype(i)) - if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) - else - itj1=ntortyp - endif - itk=itortyp(itype(k)) - itk1=itortyp(itype(k+1)) - if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) - else - itl1=ntortyp - endif +c iti=itortyp(itype(i)) +c if (j.lt.nres-1) then +c itj1=itortyp(itype(j+1)) +c else +c itj1=ntortyp +c endif +c itk=itortyp(itype(k)) +c itk1=itortyp(itype(k+1)) +c if (l.lt.nres-1) then +c itl1=itortyp(itype(l+1)) +c else +c itl1=ntortyp +c endif #ifdef MOMENT s1=dip(4,jj,i)*dip(4,kk,k) #endif @@ -9669,7 +9736,7 @@ C energy moment and not to the cluster cumulant. s2=0.5d0*scalar2(b1(1,k),auxvec(1)) call matvec2(AECA(1,1,2),b1(1,l+1),auxvec(1)) s3=0.5d0*scalar2(b1(1,j+1),auxvec(1)) - call transpose2(EE(1,1,itk),auxmat(1,1)) + call transpose2(EE(1,1,k),auxmat(1,1)) call matmat2(auxmat(1,1),AECA(1,1,1),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -9767,25 +9834,25 @@ C C 4/7/01 AL Component s1 was removed, because it pertains to the respective C energy moment and not to the cluster cumulant. cd write (2,*) 'eello_graph4: wturn6',wturn6 - iti=itortyp(itype(i)) - itj=itortyp(itype(j)) - if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) - else - itj1=ntortyp - endif - itk=itortyp(itype(k)) - if (k.lt.nres-1) then - itk1=itortyp(itype(k+1)) - else - itk1=ntortyp - endif - itl=itortyp(itype(l)) - if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) - else - itl1=ntortyp - endif +c iti=itortyp(itype(i)) +c itj=itortyp(itype(j)) +c if (j.lt.nres-1) then +c itj1=itortyp(itype(j+1)) +c else +c itj1=ntortyp +c endif +c itk=itortyp(itype(k)) +c if (k.lt.nres-1) then +c itk1=itortyp(itype(k+1)) +c else +c itk1=ntortyp +c endif +c itl=itortyp(itype(l)) +c if (l.lt.nres-1) then +c itl1=itortyp(itype(l+1)) +c else +c itl1=ntortyp +c endif cd write (2,*) 'eello6_graph4:','i',i,' j',j,' k',k,' l',l cd write (2,*) 'iti',iti,' itj',itj,' itj1',itj1,' itk',itk, cd & ' itl',itl,' itl1',itl1 @@ -10798,7 +10865,7 @@ C gradient po costhet &*VofOverlap C grad_shield_side is Cbeta sidechain gradient grad_shield_side(j,ishield_list(i),i)= - & (sh_frac_dist_grad(j)*-2.0d0 + & (sh_frac_dist_grad(j)*(-2.0d0) & +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet) & +scale_fac_dist*(cosphi_grad_long(j)) & *2.0d0/(1.0-cosphi))