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=eef70b509bab31283c5ffad560ca6176b616b8d3;hb=53a17fa4587252a39b479d21a6062bfdbe5ccb81;hp=c8bb8eba6e2786a63b75d01e8ec73db7d866b6ac;hpb=a0f6fbb2000f42d8b07a781e01f1cecc826004d1;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 c8bb8eb..eef70b5 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -85,7 +85,7 @@ C FG slaves receive the WEIGHTS array time_Bcastw=time_Bcastw+MPI_Wtime()-time00 c call chainbuild_cart endif -c print *,'Processor',myrank,' calling etotal ipot=',ipot +C print *,'Processor',myrank,' calling etotal ipot=',ipot c print *,'Processor',myrank,' nnt=',nnt,' nct=',nct #else c if (modecalc.eq.12.or.modecalc.eq.14) then @@ -98,6 +98,7 @@ c endif C C Compute the side-chain and electrostatic interaction energy C +C write(iout,*) "zaczynam liczyc energie" goto (101,102,103,104,105,106) ipot C Lennard-Jones potential. 101 call elj(evdw) @@ -117,10 +118,13 @@ C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence). goto 107 C Soft-sphere potential 106 call e_softsphere(evdw) +C write(iout,*) "skonczylem ipoty" + C C Calculate electrostatic (H-bonding) energy of the main chain. C 107 continue +C write(iout,*) "skonczylem ipoty" cmc cmc Sep-06: egb takes care of dynamic ss bonds too cmc @@ -442,7 +446,6 @@ cMS$ATTRIBUTES C :: proc_proc #endif double precision gradbufc(3,maxres),gradbufx(3,maxres), & glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres) -#endif include 'COMMON.SETUP' include 'COMMON.IOUNITS' include 'COMMON.FFIELD' @@ -715,7 +718,7 @@ c enddo do i=1,4*nres glocbuf(i)=gloc(i,icg) enddo -#define DEBUG +#undef DEBUG #ifdef DEBUG write (iout,*) "gloc_sc before reduce" do i=1,nres @@ -744,7 +747,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 +#undef DEBUG #ifdef DEBUG write (iout,*) "gloc_sc after reduce" do i=1,nres @@ -2338,6 +2341,7 @@ C endif c if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then if (i.gt. nnt+2 .and. i.lt.nct+2) then +c write(iout,*) (itype(i-2)) iti = itortyp(itype(i-2)) else iti=ntortyp+1 @@ -2391,7 +2395,11 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then do k=1,2 mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1) enddo -cd write (iout,*) 'mu ',mu(:,i-2) +cd write (iout,*) 'mu ',mu(:,i-2),i-2 +cd write (iout,*) 'b1 ',b1(:,iti1),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(:,itortyp(itype(i))),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) @@ -2718,7 +2726,7 @@ C dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), - & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4) + & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),eel_loc_ij common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 @@ -2761,6 +2769,7 @@ c call vec_and_deriv time01=MPI_Wtime() #endif call set_matrices +c write (iout,*) "after set matrices" #ifdef TIMING time_mat=time_mat+MPI_Wtime()-time01 #endif @@ -2797,6 +2806,7 @@ c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms C C Loop over i,i+2 and i,i+3 pairs of the peptide groups C +c write(iout,*) "przed turnem3 loop" do i=iturn3_start,iturn3_end if (itype(i).eq.21 .or. itype(i+1).eq.21 & .or. itype(i+2).eq.21 .or. itype(i+3).eq.21) cycle @@ -2889,7 +2899,7 @@ C------------------------------------------------------------------------------- dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), - & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4) + & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),a22,a23,a32,a33 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 @@ -3659,7 +3669,7 @@ c write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2 iti1=itortyp(itype(i+1)) iti2=itortyp(itype(i+2)) iti3=itortyp(itype(i+3)) -c write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3 +C write(iout,*) i,"iti1",iti1," iti2",iti2," iti3",iti3,itype(i+3) call transpose2(EUg(1,1,i+1),e1t(1,1)) call transpose2(Eug(1,1,i+2),e2t(1,1)) call transpose2(Eug(1,1,i+3),e3t(1,1)) @@ -4239,7 +4249,7 @@ c & "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax) else diff = vbld(i)-vbldp0 - if (energy_dec) write (iout,*) + 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 @@ -4540,7 +4550,8 @@ C logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 do i=ithet_start,ithet_end - if (itype(i-1).eq.21) cycle + if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. + &(itype(i).eq.ntyp1)) cycle dethetai=0.0d0 dephii=0.0d0 dephii1=0.0d0 @@ -4550,7 +4561,8 @@ C coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(i-2).ne.21) then +C if (i.gt.3) 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 @@ -4564,13 +4576,13 @@ C 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 enddo endif - if (i.lt.nres .and. itype(i).ne.21) then + if ((i.lt.nres).and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -4585,7 +4597,7 @@ C enddo else phii1=0.0d0 - ityp3=nthetyp+1 + ityp3=ithetyp(itype(i)) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@ -4695,6 +4707,8 @@ C & i,theta(i)*rad2deg,phii*rad2deg, & phii1*rad2deg,ethetai 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 @@ -5136,6 +5150,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 #ifdef DEBUG C @@ -5553,7 +5569,8 @@ c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end if (itype(i-2).eq.21 .or. itype(i-1).eq.21 - & .or. itype(i).eq.21) cycle + & .or. itype(i).eq.21 + & .or. itype(i-3).eq.ntyp1) cycle etors_ii=0.0D0 itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) @@ -5641,15 +5658,18 @@ 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 -c write(iout,*) "a tu??" +C write(iout,*) "a tu??" do i=iphid_start,iphid_end if (itype(i-2).eq.21 .or. itype(i-1).eq.21 - & .or. itype(i).eq.21 .or. itype(i+1).eq.21) cycle + & .or. itype(i).eq.21 .or. itype(i+1).eq.21 + & .or. itype(i-3).eq.ntyp1) cycle + etors_d_ii=0.0D0 itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) @@ -5669,6 +5689,8 @@ C Regular cosine and sine terms 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 @@ -5684,12 +5706,17 @@ C Regular cosine and sine terms 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 enddo @@ -5725,12 +5752,14 @@ c lprn=.true. c write (iout,*) "EBACK_SC_COR",itau_start,itau_end esccor=0.0D0 do i=itau_start,itau_end - 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)) 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 @@ -5754,10 +5783,14 @@ 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 -c write (iout,*) "EBACK_SC_COR",i,esccor,intertyp + if (energy_dec) write (iout,'(a6,i5,i2,0pf7.3)') + & 'esccor',i,intertyp,esccor_ii +cd write (iout,*) "tau ",i,intertyp,tauangle(intertyp,i)*RAD2DEG +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) & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')