X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fenergy_p_new_barrier.F;h=115cb57a921207c3daaf434182b7f02d7903054e;hb=582d0baf775d8a72f8fd25bac504c98b7f927d86;hp=7d54a787c6396ab292d6a05e6b236f2eef21ae69;hpb=9d5038b8ba706ed6e365dd400833fceb02b7f0a8;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 7d54a78..115cb57 100644 --- a/source/unres/src_MD/energy_p_new_barrier.F +++ b/source/unres/src_MD/energy_p_new_barrier.F @@ -263,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, @@ -822,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) @@ -1781,9 +1781,10 @@ c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 & evdwij endif - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'evdw',i,j,evdwij - + if (energy_dec) then + write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij + call flush(iout) + endif C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 fac=-expon*(e1+evdwij)*rij_shift @@ -4584,6 +4585,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) @@ -4602,6 +4605,12 @@ 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) then + write (iout,*) + & "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff, + & AKSC(1,iti),AKSC(1,iti)*diff*diff + call flush(iout) + endif 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) @@ -4894,7 +4903,7 @@ c & " ithet_end",ithet_end sinkt(k)=dsin(k*theti2) enddo C if (i.gt.3) then - if (i.gt.3 .and. itype(i-3).ne.ntyp1) then + if (i.gt.3 .and. itype(imax0(i-3,1)).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -5041,6 +5050,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 @@ -6018,7 +6029,7 @@ c c - do i=1,19 + do i=1,max_template distancek(i)=9999999.9 enddo @@ -6038,6 +6049,7 @@ 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) @@ -6058,7 +6070,18 @@ c endif enddo - min_odl=minval(distancek) + +c min_odl=minval(distancek) + do kk=1,constr_homology + if(l_homo(kk,ii)) then + min_odl=distancek(kk) + exit + endif + enddo + do kk=1,constr_homology + if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) + & min_odl=distancek(kk) + enddo c write (iout,* )"min_odl",min_odl #ifdef DEBUG write (iout,*) "ij dij",i,j,dij @@ -6071,6 +6094,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 @@ -6120,6 +6144,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 @@ -6209,7 +6234,7 @@ c write (iout,*) idihconstr_start_homo,idihconstr_end_homo do i=idihconstr_start_homo,idihconstr_end_homo kat2=0.0d0 c betai=beta(i,i+1,i+2,i+3) - betai = phi(i+3) + betai = phi(i) c write (iout,*) "betai =",betai do k=1,constr_homology dih_diff(k)=pinorm(dih(k,i)-betai) @@ -6579,12 +6604,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)) @@ -6603,6 +6630,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 @@ -6618,12 +6647,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)