From c987a84c8742686406ea6b8807c9132098dc567e Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Wed, 16 May 2018 17:26:41 +0200 Subject: [PATCH] unres src_MD-M now HOMOL energy --- ctest/prota_unres_energy_check.sh | 4 ++-- ctest/prota_unres_energy_check_mult.sh | 4 ++-- source/unres/src_MD-M/energy_p_new_barrier.F | 31 +++++++++++++++++++++----- 3 files changed, 30 insertions(+), 9 deletions(-) diff --git a/ctest/prota_unres_energy_check.sh b/ctest/prota_unres_energy_check.sh index febf018..1a8d602 100755 --- a/ctest/prota_unres_energy_check.sh +++ b/ctest/prota_unres_energy_check.sh @@ -266,8 +266,8 @@ elif [ "$1" == "remd_dfa" ]; then grep ACC remd_dfa.out_GB000 |tail -7 exchange=`grep ACC remd_dfa.out_GB000 |tail -7|awk '{a=a+$4}END{print a/NR}'` echo "average exchange = ${exchange}" - if [ `echo "a=${exchange};a<0.1"|bc -l` != "0" ]; then - echo 'ERROR average exchange smaller than 0.1' + if [ `echo "a=${exchange};a<0.05"|bc -l` != "0" ]; then + echo 'ERROR average exchange smaller than 0.05' exit 1 else exit 0 diff --git a/ctest/prota_unres_energy_check_mult.sh b/ctest/prota_unres_energy_check_mult.sh index 8831cfe..be6712e 100755 --- a/ctest/prota_unres_energy_check_mult.sh +++ b/ctest/prota_unres_energy_check_mult.sh @@ -223,8 +223,8 @@ elif [ "$1" == "Ts866_checkgrad_full" ] ; then ene=`grep ETOT $file|head -1| awk '{print $2*1.0}'` hene=`grep H_CONS $file|head -1| awk '{print $2*1.0}'` - ene_ref=3129.43 - hene_ref=1878.060 + ene_ref=3124.347 + hene_ref=1872.981 echo "ETOT= " ${ene} " reference= " ${ene_ref} if [ `echo "a=${ene}-(${ene_ref});if(0>a)a*=-1;a>0.01"|bc -l` != "0" ]; then 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 526b4ab..ea2fe9e 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -6895,7 +6895,7 @@ c c - do i=1,19 + do i=1,max_template distancek(i)=9999999.9 enddo @@ -6914,9 +6914,13 @@ c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d j = jres_homo(ii) dij=dist(i,j) c write (iout,*) "dij(",i,j,") =",dij + nexl=0 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 + if(.not.l_homo(k,ii)) then + nexl=nexl+1 + cycle + endif distance(k)=odl(k,ii)-dij c write (iout,*) "distance(",k,") =",distance(k) c @@ -6955,7 +6959,15 @@ c write (iout,* )"min_odl",min_odl write (iout,*) "distancek",(distancek(k),k=1,constr_homology) write (iout,* )"min_odl",min_odl #endif +#ifdef OLDRESTR odleg2=0.0d0 +#else + if (waga_dist.ge.0.0d0) then + odleg2=nexl + else + odleg2=0.0d0 + endif +#endif do k=1,constr_homology c Nie wiem po co to liczycie jeszcze raz! c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/ @@ -7110,8 +7122,11 @@ c if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)= c & -(6.28318-dih_diff(i,k)) c if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)= c & 6.28318+dih_diff(i,k) - +#ifdef OLD_DIHED kat3=-0.5d0*dih_diff(k)**2*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument +#else + kat3=(dcos(dih_diff(k))-1)*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument +#endif c kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i) gdih(k)=dexp(kat3) kat2=kat2+gdih(k) @@ -7138,7 +7153,11 @@ c ---------------------------------------------------------------------- sum_gdih=kat2 sum_sgdih=0.0d0 do k=1,constr_homology +#ifdef OLD_DIHED sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i) ! waga_angle rmvd +#else + sgdih=-gdih(k)*dsin(dih_diff(k))*sigma_dih(k,i) ! waga_angle rmvd +#endif c sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle sum_sgdih=sum_sgdih+sgdih enddo @@ -7220,7 +7239,7 @@ c utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument c utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta? gtheta(k)=dexp(utheta_i) ! + min_utheta_i? - gutheta_i=gutheta_i+dexp(utheta_i) ! Sum of Gaussians (pk) + gutheta_i=gutheta_i+gtheta(k) ! Sum of Gaussians (pk) c Gradient for single Gaussian restraint in subr Econstr_back c dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*dtheta_i/(ii-1) c @@ -7292,7 +7311,9 @@ c c usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d? c uscdiffk(k)=usc_diff(i) guscdiff2(k)=dexp(usc_diff_i) ! without min_scdiff - guscdiff(i)=guscdiff(i)+dexp(usc_diff_i) !Sum of Gaussians (pk) +c write(iout,*) "i",i," k",k," sigma_d",sigma_d(k,i), +c & " guscdiff2",guscdiff2(k) + guscdiff(i)=guscdiff(i)+guscdiff2(k) !Sum of Gaussians (pk) c write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j), c & xxref(j),yyref(j),zzref(j) enddo -- 1.7.9.5