X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fenergy_p_new_barrier.F;h=f40a9930ef73ed1659ad20ecbef787f59fcda52f;hb=37982ff2ef784ba689f3f1c1cbf3337d9daa0cd9;hp=002f211b43fbd6eb4c9d219f8d0aefcf8a56ee09;hpb=aa4540154949fdfc9ab63442eb6fea09d6c589ad;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 002f211..f40a993 100644 --- a/source/unres/src_MD/energy_p_new_barrier.F +++ b/source/unres/src_MD/energy_p_new_barrier.F @@ -99,6 +99,12 @@ c if (modecalc.eq.12.or.modecalc.eq.14) then c call int_from_cart1(.false.) c endif #endif +#ifndef DFA + edfadis=0.0d0 + edfator=0.0d0 + edfanei=0.0d0 + edfabet=0.0d0 +#endif #ifdef TIMING #ifdef MPI time00=MPI_Wtime() @@ -132,6 +138,7 @@ C C Calculate electrostatic (H-bonding) energy of the main chain. C 107 continue +#ifdef DFA C BARTEK for dfa test! if (wdfa_dist.gt.0) then call edfad(edfadis) @@ -156,6 +163,7 @@ c print*, 'edfan is finished!', edfanei else edfabet=0 endif +#endif c print*, 'edfab is finished!', edfabet cmc cmc Sep-06: egb takes care of dynamic ss bonds too @@ -255,7 +263,7 @@ cd print *,'nterm=',nterm edihcnstr=0 endif - if (constr_homology.ge.1) then + if (constr_homology.ge.1.and.waga_homology(iset).ne.0d0) then call e_modeller(ehomology_constr) c print *,'iset=',iset,'me=',me,ehomology_constr, c & 'Processor',fg_rank,' CG group',kolor, @@ -814,7 +822,7 @@ c enddo #endif enddo enddo - if (constr_homology.gt.0) then + if (constr_homology.gt.0.and.waga_homology(iset).ne.0d0) then do i=1,nct do j=1,3 gradc(j,i,icg)=gradc(j,i,icg)+duscdiff(j,i) @@ -4576,6 +4584,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) @@ -4594,6 +4604,9 @@ 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) write (iout,*) + & "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff, + & AKSC(1,iti),AKSC(1,iti)*diff*diff 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) @@ -5033,6 +5046,8 @@ c lprn1=.true. & 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)=gloc(nphi+i-2,icg)+wang*dethetai @@ -6010,7 +6025,7 @@ c c - do i=1,19 + do i=1,max_template distancek(i)=9999999.9 enddo @@ -6030,6 +6045,8 @@ c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d dij=dist(i,j) c write (iout,*) "dij(",i,j,") =",dij do k=1,constr_homology +c write(iout,*) ii,k,i,j,l_homo(k,ii),dij,odl(k,ii) + if(.not.l_homo(k,ii)) cycle distance(k)=odl(k,ii)-dij c write (iout,*) "distance(",k,") =",distance(k) c @@ -6062,6 +6079,7 @@ c write (iout,* )"min_odl",min_odl c Nie wiem po co to liczycie jeszcze raz! c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/ c & (2*(sigma_odl(i,j,k))**2)) + if(.not.l_homo(k,ii)) cycle if (waga_dist.ge.0.0d0) then c c For Gaussian-type Urestr @@ -6111,6 +6129,7 @@ c godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2)) c & *waga_dist)+min_odl c sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist c + if(.not.l_homo(k,ii)) cycle if (waga_dist.ge.0.0d0) then c For Gaussian-type Urestr c @@ -6335,9 +6354,6 @@ c c sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i)*waga_theta ! right functional form? sum_sgtheta=sum_sgtheta+sgtheta ! cum variable enddo -c grad_theta3=sum_sgtheta/sum_gtheta 1/*theta(i)? s. line below -c grad_theta3=sum_sgtheta/sum_gtheta -c c Final value of gradient using same var as in Econstr_back gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg) & +sum_sgtheta/sum_gtheta*waga_theta @@ -6573,12 +6589,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)) @@ -6597,6 +6615,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 @@ -6612,12 +6632,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)