X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fenergy_p_new_barrier.F;h=e5949761a5d67a6e297fa446c0410cfba9cb681a;hb=132235b8d133036519f7e98dd9dde507de4958a8;hp=d07d13553c2629b5e9ffc61886eef2b70228487c;hpb=989b55752ffaadfa073c68b17130a1755d304d6b;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 d07d135..e594976 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 @@ -543,6 +551,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.TIME1' include 'COMMON.MAXGRAD' include 'COMMON.SCCOR' + include 'COMMON.MD' #ifdef TIMING #ifdef MPI time01=MPI_Wtime() @@ -813,6 +822,14 @@ c enddo #endif enddo enddo + if (constr_homology.gt.0) then + do i=1,nct + do j=1,3 + gradc(j,i,icg)=gradc(j,i,icg)+duscdiff(j,i) + gradx(j,i,icg)=gradx(j,i,icg)+duscdiffx(j,i) + enddo + enddo + endif #ifdef DEBUG write (iout,*) "gloc before adding corr" do i=1,4*nres @@ -4862,6 +4879,8 @@ C & sinph1ph2(maxdouble,maxdouble) logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 +c write (iout,*) "EBEND ithet_start",ithet_start, +c & " ithet_end",ithet_end do i=ithet_start,ithet_end if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. &(itype(i).eq.ntyp1)) cycle @@ -6019,6 +6038,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 @@ -6051,6 +6072,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 @@ -6100,6 +6122,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 @@ -6324,12 +6347,12 @@ 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 - dutheta(i-2)=sum_sgtheta/sum_gtheta*waga_theta + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg) + & +sum_sgtheta/sum_gtheta*waga_theta & *waga_homology(iset) +c dutheta(i-2)=sum_sgtheta/sum_gtheta*waga_theta +c & *waga_homology(iset) c dutheta(i)=sum_sgtheta/sum_gtheta c c Uconst_back=Uconst_back+waga_theta*utheta(i) ! waga_theta added as weight @@ -6525,9 +6548,12 @@ c & waga_theta*Eval+waga_d*Erot)*waga_homology(iset) c write (iout,*) "ehomology_constr=",ehomology_constr endif -c write (iout,*) "odleg",odleg," kat",kat," Uconst_back",Uconst_back -c write (iout,*) "ehomology_constr",ehomology_constr -c ehomology_constr=odleg+kat+Uconst_back +#ifdef DEBUG + write (iout,*) "odleg",waga_dist,odleg," kat",waga_angle,kat, + & "Eval",waga_theta,eval, + & "Erot",waga_d,Erot + write (iout,*) "ehomology_constr",ehomology_constr +#endif return c c FP 01/15 end @@ -6556,12 +6582,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)) @@ -6580,6 +6608,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 @@ -6595,12 +6625,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)