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=f95ca0c840122660eb1463b376658b093eee9571;hb=024c09f85d96670bfddf91ec2305f8a5e2532adc;hp=b3c4c22ff9bc23562f5512887d31cc3f9407927a;hpb=3c5e1a8829ff84c0725642015445f4debd4ecc0a;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 b3c4c22..f95ca0c 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -55,6 +55,7 @@ C FG slaves as WEIGHTS array. weights_(17)=wbond weights_(18)=scal14 weights_(21)=wsccor + weights_(25)=wsaxs C FG Master broadcasts the WEIGHTS_ array call MPI_Bcast(weights_(1),n_ene, & MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR) @@ -81,6 +82,7 @@ C FG slaves receive the WEIGHTS array wbond=weights(17) scal14=weights(18) wsccor=weights(21) + wsaxs=weights(25) endif time_Bcast=time_Bcast+MPI_Wtime()-time00 time_Bcastw=time_Bcastw+MPI_Wtime()-time00 @@ -266,6 +268,16 @@ cd &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6 call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1) cd write (iout,*) "multibody_hb ecorr",ecorr endif +c write (iout,*) "nsaxs",nsaxs," saxs_mode",saxs_mode + if (nsaxs.gt.0 .and. saxs_mode.eq.0) then + call e_saxs(Esaxs_constr) +c write (iout,*) "From Esaxs: Esaxs_constr",Esaxs_constr + else if (nsaxs.gt.0 .and. saxs_mode.gt.0) then + call e_saxsC(Esaxs_constr) +c write (iout,*) "From EsaxsC: Esaxs_constr",Esaxs_constr + else + Esaxs_constr = 0.0d0 + endif c print *,"Processor",myrank," computed Ucorr" C C If performing constraint dynamics, call the constraint energy @@ -334,6 +346,7 @@ C energia(22)=eliptran energia(23)=Eafmforce energia(24)=ehomology_constr + energia(25)=Esaxs_constr c Here are the energies showed per procesor if the are more processors c per molecule then we sum it up in sum_energy subroutine c print *," Processor",myrank," calls SUM_ENERGY" @@ -428,6 +441,9 @@ cMS$ATTRIBUTES C :: proc_proc eliptran=energia(22) Eafmforce=energia(23) ehomology_constr=energia(24) + esaxs_constr=energia(25) +c write (iout,*) "sum_energy esaxs_constr",esaxs_constr, +c & " wsaxs",wsaxs #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc @@ -435,7 +451,7 @@ cMS$ATTRIBUTES C :: proc_proc & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+Uconst+wsccor*esccor+ehomology_constr - & +wliptran*eliptran+Eafmforce + & +wsaxs*esaxs_constr+wliptran*eliptran+Eafmforce #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc @@ -443,7 +459,7 @@ cMS$ATTRIBUTES C :: proc_proc & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+Uconst+wsccor*esccor+ehomology_constr - & +wliptran*eliptran + & +wsaxs*esaxs_constr+wliptran*eliptran & +Eafmforce #endif energia(0)=etot @@ -502,12 +518,20 @@ cMS$ATTRIBUTES C :: proc_proc #endif #ifdef DEBUG write (iout,*) "sum_gradient gvdwc, gvdwx" - do i=1,nres - write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') + do i=0,nres + write (iout,'(i3,3e15.5,5x,3e15.5)') & i,(gvdwx(j,i),j=1,3),(gvdwc(j,i),j=1,3) enddo call flush(iout) #endif +#ifdef DEBUG + write (iout,*) "sum_gradient gsaxsc, gsaxsx" + do i=0,nres + write (iout,'(i3,3e15.5,5x,3e15.5)') + & i,(gsaxsc(j,i),j=1,3),(gsaxsx(j,i),j=1,3) + enddo + call flush(iout) +#endif #ifdef MPI C FG slaves call the following matching MPI_Bcast in ERGASTULUM if (nfgtasks.gt.1 .and. fg_rank.eq.0) @@ -547,7 +571,8 @@ c enddo & wcorr5*gradcorr5_long(j,i)+ & wcorr6*gradcorr6_long(j,i)+ & wturn6*gcorr6_turn_long(j,i)+ - & wstrain*ghpbc(j,i) + & wstrain*ghpbc(j,i)+ + & wsaxs*gsaxsc(j,i) & +wliptran*gliptranc(j,i) & +gradafm(j,i) @@ -565,7 +590,8 @@ c enddo & wcorr5*gradcorr5_long(j,i)+ & wcorr6*gradcorr6_long(j,i)+ & wturn6*gcorr6_turn_long(j,i)+ - & wstrain*ghpbc(j,i) + & wstrain*ghpbc(j,i)+ + & wsaxs*gsaxsc(j,i) & +wliptran*gliptranc(j,i) & +gradafm(j,i) @@ -625,7 +651,8 @@ c enddo do j=1,3 gradbufc(j,nres-1)=gradbufc_sum(j,nres) enddo - do i=nres-2,-1,-1 + do i=nres-2,0,-1 +c do i=nres-2,nnt,-1 do j=1,3 gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1) enddo @@ -633,7 +660,7 @@ c enddo #ifdef DEBUG write (iout,*) "gradbufc after summing" do i=1,nres - write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3) + write (iout,'(i3,3e15.5)') i,(gradbufc(j,i),j=1,3) enddo call flush(iout) #endif @@ -641,8 +668,8 @@ c enddo #endif #ifdef DEBUG write (iout,*) "gradbufc" - do i=1,nres - write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3) + do i=0,nres + write (iout,'(i3,3e15.5)') i,(gradbufc(j,i),j=1,3) enddo call flush(iout) #endif @@ -655,7 +682,8 @@ c enddo do j=1,3 gradbufc(j,nres-1)=gradbufc_sum(j,nres) enddo - do i=nres-2,-1,-1 + do i=nres-2,0,-1 +c do i=nres-2,nnt,-1 do j=1,3 gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1) enddo @@ -672,8 +700,8 @@ c enddo c enddo #ifdef DEBUG write (iout,*) "gradbufc after summing" - do i=1,nres - write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3) + do i=0,nres + write (iout,'(i3,3e15.5)') i,(gradbufc(j,i),j=1,3) enddo call flush(iout) #endif @@ -731,8 +759,9 @@ c enddo #endif gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ - & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ - & wsccor*gsccorx(j,i) + & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i) + & +wsaxs*gsaxsx(j,i) + & +wsccor*gsccorx(j,i) & +wscloc*gsclocx(j,i) & +wliptran*gliptranx(j,i) enddo @@ -929,7 +958,7 @@ c #ifdef DEBUG write (iout,*) "gradc gradx gloc" do i=1,nres - write (iout,'(i5,3f10.5,5x,3f10.5,5x,f10.5)') + write (iout,'(i5,3e15.5,5x,3e15.5,5x,e15.5)') & i,(gradc(j,i,icg),j=1,3),(gradx(j,i,icg),j=1,3),gloc(i,icg) enddo #endif @@ -1032,6 +1061,7 @@ C------------------------------------------------------------------------ Uconst=energia(20) esccor=energia(21) ehomology_constr=energia(24) + esaxs_constr=energia(25) eliptran=energia(22) Eafmforce=energia(23) #ifdef SPLITELE @@ -1041,7 +1071,7 @@ C------------------------------------------------------------------------ & ecorr,wcorr, & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor, - & edihcnstr,ehomology_constr, ebr*nss, + & edihcnstr,ehomology_constr,esaxs_constr*wsaxs, ebr*nss, & Uconst,eliptran,wliptran,Eafmforce,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ @@ -1065,6 +1095,7 @@ C------------------------------------------------------------------------ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/ + & 'E_SAXS=',1pE16.6,' (SAXS restraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST= ',1pE16.6,' (Constraint energy)'/ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ @@ -1078,7 +1109,7 @@ C------------------------------------------------------------------------ & ecorr,wcorr, & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr, - & ehomology_constr,ebr*nss,Uconst, + & ehomology_constr,esaxs_constr*wsaxs,ebr*nss,Uconst, & eliptran,wliptran,Eafmforc, & etot 10 format (/'Virtual-chain energies:'// @@ -1102,6 +1133,7 @@ C------------------------------------------------------------------------ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/ + & 'E_SAXS=',1pE16.6,' (SAXS restraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST=',1pE16.6,' (Constraint energy)'/ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ @@ -2860,7 +2892,8 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then cd write (iout,*) '*******i',i,' iti1',iti cd write (iout,*) 'b1',b1(:,iti) cd write (iout,*) 'b2',b2(:,iti) -cd write (iout,*) 'Ug',Ug(:,:,i-2) +cd write (iout,*) "phi(",i,")=",phi(i)," sin1",sin1," cos1",cos1 +cd write (iout,*) 'Ug',Ug(:,:,i-2) c if (i .gt. iatel_s+2) then if (i .gt. nnt+2) then call matvec2(Ug(1,1,i-2),b2(1,i-2),Ub2(1,i-2)) @@ -2915,7 +2948,11 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1) enddo C write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2) -c write (iout,*) 'mu ',mu(:,i-2),i-2 +cd write (iout,*) 'mu ',mu(:,i-2),i-2 +cd write (iout,*) 'b1 ',b1(:,i-1),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(:,i-2),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) @@ -3324,21 +3361,21 @@ C Loop over i,i+2 and i,i+3 pairs of the peptide groups C C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition do i=iturn3_start,iturn3_end - if (i.le.1) cycle +CAna if (i.le.1) cycle C write(iout,*) "tu jest i",i if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((i+4).gt.nres) - & .or.((i-1).le.0) +CAna & .or.((i+4).gt.nres) +CAna & .or.((i-1).le.0) C end of changes by Ana & .or. itype(i+2).eq.ntyp1 & .or. itype(i+3).eq.ntyp1) cycle - if(i.gt.1)then - if(itype(i-1).eq.ntyp1)cycle - end if - if(i.LT.nres-3)then - if (itype(i+4).eq.ntyp1) cycle - end if +CAna if(i.gt.1)then +CAna if(itype(i-1).eq.ntyp1)cycle +CAna end if +CAna if(i.LT.nres-3)then +CAna if (itype(i+4).eq.ntyp1) cycle +CAna end if dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -3360,17 +3397,17 @@ C end of changes by Ana num_cont_hb(i)=num_conti enddo do i=iturn4_start,iturn4_end - if (i.le.1) cycle +cAna if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((i+5).gt.nres) - & .or.((i-1).le.0) +cAna & .or.((i+5).gt.nres) +cAna & .or.((i-1).le.0) C end of changes suggested by Ana & .or. itype(i+3).eq.ntyp1 & .or. itype(i+4).eq.ntyp1 - & .or. itype(i+5).eq.ntyp1 - & .or. itype(i).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +cAna & .or. itype(i+5).eq.ntyp1 +cAna & .or. itype(i).eq.ntyp1 +cAna & .or. itype(i-1).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -3428,14 +3465,14 @@ c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c do i=iatel_s,iatel_e - if (i.le.1) cycle +cAna if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((i+2).gt.nres) - & .or.((i-1).le.0) +cAna & .or.((i+2).gt.nres) +cAna & .or.((i-1).le.0) C end of changes by Ana - & .or. itype(i+2).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +cAna & .or. itype(i+2).eq.ntyp1 +cAna & .or. itype(i-1).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -3486,14 +3523,14 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) do j=ielstart(i),ielend(i) C write (iout,*) i,j - if (j.le.1) cycle +cAna if (j.le.1) cycle if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((j+2).gt.nres) - & .or.((j-1).le.0) +cAna & .or.((j+2).gt.nres) +cAna & .or.((j-1).le.0) C end of changes by Ana - & .or.itype(j+2).eq.ntyp1 - & .or.itype(j-1).eq.ntyp1 +cAna & .or.itype(j+2).eq.ntyp1 +cAna & .or.itype(j-1).eq.ntyp1 &) cycle call eelecij(i,j,ees,evdw1,eel_loc) enddo ! j @@ -4388,6 +4425,7 @@ C Derivatives in gamma(i+1) gel_loc_turn3(i+1)=gel_loc_turn3(i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) C Cartesian derivatives +!DIR$ UNROLL(0) do l=1,3 c ghalf1=0.5d0*agg(l,1) c ghalf2=0.5d0*agg(l,2) @@ -5684,7 +5722,7 @@ C logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 do i=ithet_start,ithet_end - if (i.eq.2) cycle +c if (i.eq.2) cycle c print *,i,itype(i-1),itype(i),itype(i-2) if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1) & .or.(itype(i).eq.ntyp1)) cycle @@ -5701,7 +5739,7 @@ C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(i-3).ne.ntyp1) 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 @@ -5716,7 +5754,7 @@ C propagation of chirality for glycine type 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 @@ -5737,7 +5775,7 @@ C propagation of chirality for glycine type enddo else phii1=0.0d0 - ityp3=nthetyp+1 + ityp3=ithetyp(itype(i)) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@ -6295,6 +6333,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,it,itype(i) c & ,zz,xx,yy c#define DEBUG @@ -7060,11 +7100,12 @@ 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) -c write (iout,*) "dih_diff(",k,") =",dih_diff(k) +cd write (iout,'(a8,2i4,2f15.8)') "dih_diff",i,k,dih_diff(k) +cd & ,sigma_dih(k,i) 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)= @@ -7107,7 +7148,7 @@ c grad_dih3=sum_sgdih/sum_gdih c write(iout,*)i,k,gdih,sgdih,beta(i+1,i+2,i+3,i+4),grad_dih3 ccc write(iout,747) "GRAD_KAT_1", i, nphi, icg, grad_dih3, ccc & gloc(nphi+i-3,icg) - gloc(i,icg)=gloc(i,icg)+grad_dih3 + gloc(i-3,icg)=gloc(i-3,icg)+grad_dih3 c if (i.eq.25) then c write(iout,*) "i",i,"icg",icg,"gloc(",i,icg,")",gloc(i,icg) c endif @@ -7172,6 +7213,9 @@ c c dtheta_i=theta(j)-thetaref(j,iref) c dtheta_i=thetaref(k,i)-theta(i) ! original form without indexing theta_diff(k)=thetatpl(k,i)-theta(i) +cd write (iout,'(a8,2i4,2f15.8)') "theta_diff",i,k,theta_diff(k) +cd & ,sigma_theta(k,i) + 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? @@ -7242,7 +7286,7 @@ c Original sign inverted for calc of gradient dyy=-yytpl(k,i)+yytab(i) ! ibid y dzz=-zztpl(k,i)+zztab(i) ! ibid z c write(iout,*) "dxx, dyy, dzz" -c write(iout,'(2i5,3f8.2)') k,i,dxx,dyy,dzz +cd write(iout,'(2i5,4f8.2)') k,i,dxx,dyy,dzz,sigma_d(k,i) c usc_diff_i=-0.5d0*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d rmvd from Gaussian argument c usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d? @@ -7546,12 +7590,12 @@ c write (iout,*) "EBACK_SC_COR",itau_start,itau_end esccor=0.0D0 do i=itau_start,itau_end if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle - esccor_ii=0.0D0 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 @@ -7575,9 +7619,12 @@ 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 + if (energy_dec) write (iout,'(a6,i5,i2,0pf7.3)') + & 'esccor',i,intertyp,esccor_ii 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) @@ -11002,3 +11049,403 @@ C print *,'AFM',Eafmforce,totTafm*velAFMconst,dist return end +c---------------------------------------------------------------------------- + double precision function sscale2(r,r_cut,r0,rlamb) + implicit none + double precision r,gamm,r_cut,r0,rlamb,rr + rr = dabs(r-r0) +c write (2,*) "r",r," r_cut",r_cut," r0",r0," rlamb",rlamb +c write (2,*) "rr",rr + if(rr.lt.r_cut-rlamb) then + sscale2=1.0d0 + else if(rr.le.r_cut.and.rr.ge.r_cut-rlamb) then + gamm=(rr-(r_cut-rlamb))/rlamb + sscale2=1.0d0+gamm*gamm*(2*gamm-3.0d0) + else + sscale2=0d0 + endif + return + end +C----------------------------------------------------------------------- + double precision function sscalgrad2(r,r_cut,r0,rlamb) + implicit none + double precision r,gamm,r_cut,r0,rlamb,rr + rr = dabs(r-r0) + if(rr.lt.r_cut-rlamb) then + sscalgrad2=0.0d0 + else if(rr.le.r_cut.and.rr.ge.r_cut-rlamb) then + gamm=(rr-(r_cut-rlamb))/rlamb + if (r.ge.r0) then + sscalgrad2=gamm*(6*gamm-6.0d0)/rlamb + else + sscalgrad2=-gamm*(6*gamm-6.0d0)/rlamb + endif + else + sscalgrad2=0.0d0 + endif + return + end +c---------------------------------------------------------------------------- + subroutine e_saxs(Esaxs_constr) + implicit none + include 'DIMENSIONS' +#ifdef MPI + include "mpif.h" + include "COMMON.SETUP" + integer IERR +#endif + include 'COMMON.SBRIDGE' + include 'COMMON.CHAIN' + include 'COMMON.GEO' + include 'COMMON.DERIV' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.MD' + include 'COMMON.CONTROL' + include 'COMMON.NAMES' + include 'COMMON.TIME1' + include 'COMMON.FFIELD' +c + double precision Esaxs_constr + integer i,iint,j,k,l + double precision PgradC(maxSAXS,3,maxres), + & PgradX(maxSAXS,3,maxres),Pcalc(maxSAXS) +#ifdef MPI + double precision PgradC_(maxSAXS,3,maxres), + & PgradX_(maxSAXS,3,maxres),Pcalc_(maxSAXS) +#endif + double precision dk,dijCACA,dijCASC,dijSCCA,dijSCSC, + & sigma2CACA,sigma2CASC,sigma2SCCA,sigma2SCSC,expCACA,expCASC, + & expSCCA,expSCSC,CASCgrad,SCCAgrad,SCSCgrad,aux,auxC,auxC1, + & auxX,auxX1,CACAgrad,Cnorm + double precision sss2,ssgrad2,rrr,sscalgrad2,sscale2 + double precision dist + external dist +c SAXS restraint penalty function +#ifdef DEBUG + write(iout,*) "------- SAXS penalty function start -------" + write (iout,*) "nsaxs",nsaxs + write (iout,*) "Esaxs: iatsc_s",iatsc_s," iatsc_e",iatsc_e + write (iout,*) "Psaxs" + do i=1,nsaxs + write (iout,'(i5,e15.5)') i, Psaxs(i) + enddo +#endif + Esaxs_constr = 0.0d0 + do k=1,nsaxs + Pcalc(k)=0.0d0 + do j=1,nres + do l=1,3 + PgradC(k,l,j)=0.0d0 + PgradX(k,l,j)=0.0d0 + enddo + enddo + enddo + do i=iatsc_s,iatsc_e + if (itype(i).eq.ntyp1) cycle + do iint=1,nint_gr(i) + do j=istart(i,iint),iend(i,iint) + if (itype(j).eq.ntyp1) cycle +#ifdef ALLSAXS + dijCACA=dist(i,j) + dijCASC=dist(i,j+nres) + dijSCCA=dist(i+nres,j) + dijSCSC=dist(i+nres,j+nres) + sigma2CACA=2.0d0/(pstok**2) + sigma2CASC=4.0d0/(pstok**2+restok(itype(j))**2) + sigma2SCCA=4.0d0/(pstok**2+restok(itype(i))**2) + sigma2SCSC=4.0d0/(restok(itype(j))**2+restok(itype(i))**2) + do k=1,nsaxs + dk = distsaxs(k) + expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2) + if (itype(j).ne.10) then + expCASC = dexp(-0.5d0*sigma2CASC*(dijCASC-dk)**2) + else + endif + expCASC = 0.0d0 + if (itype(i).ne.10) then + expSCCA = dexp(-0.5d0*sigma2SCCA*(dijSCCA-dk)**2) + else + expSCCA = 0.0d0 + endif + if (itype(i).ne.10 .and. itype(j).ne.10) then + expSCSC = dexp(-0.5d0*sigma2SCSC*(dijSCSC-dk)**2) + else + expSCSC = 0.0d0 + endif + Pcalc(k) = Pcalc(k)+expCACA+expCASC+expSCCA+expSCSC +#ifdef DEBUG + write(iout,*) "i j k Pcalc",i,j,Pcalc(k) +#endif + CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA + CASCgrad = sigma2CASC*(dijCASC-dk)*expCASC + SCCAgrad = sigma2SCCA*(dijSCCA-dk)*expSCCA + SCSCgrad = sigma2SCSC*(dijSCSC-dk)*expSCSC + do l=1,3 +c CA CA + aux = CACAgrad*(C(l,j)-C(l,i))/dijCACA + PgradC(k,l,i) = PgradC(k,l,i)-aux + PgradC(k,l,j) = PgradC(k,l,j)+aux +c CA SC + if (itype(j).ne.10) then + aux = CASCgrad*(C(l,j+nres)-C(l,i))/dijCASC + PgradC(k,l,i) = PgradC(k,l,i)-aux + PgradC(k,l,j) = PgradC(k,l,j)+aux + PgradX(k,l,j) = PgradX(k,l,j)+aux + endif +c SC CA + if (itype(i).ne.10) then + aux = SCCAgrad*(C(l,j)-C(l,i+nres))/dijSCCA + PgradX(k,l,i) = PgradX(k,l,i)-aux + PgradC(k,l,i) = PgradC(k,l,i)-aux + PgradC(k,l,j) = PgradC(k,l,j)+aux + endif +c SC SC + if (itype(i).ne.10 .and. itype(j).ne.10) then + aux = SCSCgrad*(C(l,j+nres)-C(l,i+nres))/dijSCSC + PgradC(k,l,i) = PgradC(k,l,i)-aux + PgradC(k,l,j) = PgradC(k,l,j)+aux + PgradX(k,l,i) = PgradX(k,l,i)-aux + PgradX(k,l,j) = PgradX(k,l,j)+aux + endif + enddo ! l + enddo ! k +#else + dijCACA=dist(i,j) + sigma2CACA=scal_rad**2*0.25d0/ + & (restok(itype(j))**2+restok(itype(i))**2) + rrr = saxs_cutoff*2.0d0/dsqrt(sigma2CACA) + do k=1,nsaxs + dk = distsaxs(k) +c write (2,*) "ijk",i,j,k + sss2 = sscale2(dijCACA,rrr,dk,0.3d0) + if (sss2.ne.0.0d0) then + ssgrad2 = sscalgrad2(dijCACA,rrr,dk,0.3d0) + if (energy_dec) write(iout,'(a4,3i5,5f10.4)') + & 'saxs',i,j,k,dijCACA,rrr,dk,sss2,ssgrad2 + expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)*sss2 + Pcalc(k) = Pcalc(k)+expCACA +#ifdef DEBUG + write(iout,*) "i j k Pcalc",i,j,Pcalc(k) +#endif + CACAgrad = -sigma2CACA*(dijCACA-dk)*expCACA+ + & ssgrad2*expCACA/sss2 + do l=1,3 +c CA CA + aux = CACAgrad*(C(l,j)-C(l,i))/dijCACA + PgradC(k,l,i) = PgradC(k,l,i)+aux + PgradC(k,l,j) = PgradC(k,l,j)-aux + enddo ! l + endif ! sss + enddo ! k +#endif + enddo ! j + enddo ! iint + enddo ! i +#ifdef MPI + if (nfgtasks.gt.1) then + call MPI_Reduce(Pcalc(1),Pcalc_(1),nsaxs,MPI_DOUBLE_PRECISION, + & MPI_SUM,king,FG_COMM,IERR) + if (fg_rank.eq.king) then + do k=1,nsaxs + Pcalc(k) = Pcalc_(k) + enddo + endif + call MPI_Reduce(PgradC(k,1,1),PgradC_(k,1,1),3*maxsaxs*nres, + & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) + if (fg_rank.eq.king) then + do i=1,nres + do l=1,3 + do k=1,nsaxs + PgradC(k,l,i) = PgradC_(k,l,i) + enddo + enddo + enddo + endif +#ifdef ALLSAXS + call MPI_Reduce(PgradX(k,1,1),PgradX_(k,1,1),3*maxsaxs*nres, + & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) + if (fg_rank.eq.king) then + do i=1,nres + do l=1,3 + do k=1,nsaxs + PgradX(k,l,i) = PgradX_(k,l,i) + enddo + enddo + enddo + endif +#endif + endif +#endif +#ifdef MPI + if (fg_rank.eq.king) then +#endif + Cnorm = 0.0d0 + do k=1,nsaxs + Cnorm = Cnorm + Pcalc(k) + enddo + Esaxs_constr = dlog(Cnorm) + do k=1,nsaxs + if (Pcalc(k).gt.0.0d0) + & Esaxs_constr = Esaxs_constr - Psaxs(k)*dlog(Pcalc(k)) +#ifdef DEBUG + write (iout,*) "k",k," Esaxs_constr",Esaxs_constr +#endif + enddo +#ifdef DEBUG + write (iout,*) "Cnorm",Cnorm," Esaxs_constr",Esaxs_constr +#endif + do i=nnt,nct + do l=1,3 + auxC=0.0d0 + auxC1=0.0d0 + auxX=0.0d0 + auxX1=0.d0 + do k=1,nsaxs + if (Pcalc(k).gt.0) + & auxC = auxC +Psaxs(k)*PgradC(k,l,i)/Pcalc(k) + auxC1 = auxC1+PgradC(k,l,i) +#ifdef ALLSAXS + auxX = auxX +Psaxs(k)*PgradX(k,l,i)/Pcalc(k) + auxX1 = auxX1+PgradX(k,l,i) +#endif + enddo + gsaxsC(l,i) = auxC - auxC1/Cnorm +#ifdef ALLSAXS + gsaxsX(l,i) = auxX - auxX1/Cnorm +#endif +c write (iout,*) "l i",l,i," gradC",wsaxs*(auxC - auxC1/Cnorm), +c * " gradX",wsaxs*(auxX - auxX1/Cnorm) + enddo + enddo +#ifdef MPI + endif +#endif + return + end +c---------------------------------------------------------------------------- + subroutine e_saxsC(Esaxs_constr) + implicit none + include 'DIMENSIONS' +#ifdef MPI + include "mpif.h" + include "COMMON.SETUP" + integer IERR +#endif + include 'COMMON.SBRIDGE' + include 'COMMON.CHAIN' + include 'COMMON.GEO' + include 'COMMON.DERIV' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.MD' + include 'COMMON.CONTROL' + include 'COMMON.NAMES' + include 'COMMON.TIME1' + include 'COMMON.FFIELD' +c + double precision Esaxs_constr + integer i,iint,j,k,l + double precision PgradC(3,maxres),PgradX(3,maxres),Pcalc,logPtot +#ifdef MPI + double precision gsaxsc_(3,maxres),gsaxsx_(3,maxres),logPtot_ +#endif + double precision dk,dijCASPH,dijSCSPH, + & sigma2CA,sigma2SC,expCASPH,expSCSPH, + & CASPHgrad,SCSPHgrad,aux,auxC,auxC1, + & auxX,auxX1,Cnorm +c SAXS restraint penalty function +#ifdef DEBUG + write(iout,*) "------- SAXS penalty function start -------" + write (iout,*) "nsaxs",nsaxs + + do i=nnt,nct + print *,MyRank,"C",i,(C(j,i),j=1,3) + enddo + do i=nnt,nct + print *,MyRank,"CSaxs",i,(Csaxs(j,i),j=1,3) + enddo +#endif + Esaxs_constr = 0.0d0 + logPtot=0.0d0 + do j=isaxs_start,isaxs_end + Pcalc=0.0d0 + do i=1,nres + do l=1,3 + PgradC(l,i)=0.0d0 + PgradX(l,i)=0.0d0 + enddo + enddo + do i=nnt,nct + if (itype(i).eq.ntyp1) cycle + dijCASPH=0.0d0 + dijSCSPH=0.0d0 + do l=1,3 + dijCASPH=dijCASPH+(C(l,i)-Csaxs(l,j))**2 + enddo + if (itype(i).ne.10) then + do l=1,3 + dijSCSPH=dijSCSPH+(C(l,i+nres)-Csaxs(l,j))**2 + enddo + endif + sigma2CA=2.0d0/pstok**2 + sigma2SC=4.0d0/restok(itype(i))**2 + expCASPH = dexp(-0.5d0*sigma2CA*dijCASPH) + expSCSPH = dexp(-0.5d0*sigma2SC*dijSCSPH) + Pcalc = Pcalc+expCASPH+expSCSPH +#ifdef DEBUG + write(*,*) "processor i j Pcalc", + & MyRank,i,j,dijCASPH,dijSCSPH, Pcalc +#endif + CASPHgrad = sigma2CA*expCASPH + SCSPHgrad = sigma2SC*expSCSPH + do l=1,3 + aux = (C(l,i+nres)-Csaxs(l,j))*SCSPHgrad + PgradX(l,i) = PgradX(l,i) + aux + PgradC(l,i) = PgradC(l,i)+(C(l,i)-Csaxs(l,j))*CASPHgrad+aux + enddo ! l + enddo ! i + do i=nnt,nct + do l=1,3 + gsaxsc(l,i)=gsaxsc(l,i)+PgradC(l,i)/Pcalc + gsaxsx(l,i)=gsaxsx(l,i)+PgradX(l,i)/Pcalc + enddo + enddo + logPtot = logPtot - dlog(Pcalc) +c print *,"me",me,MyRank," j",j," logPcalc",-dlog(Pcalc), +c & " logPtot",logPtot + enddo ! j +#ifdef MPI + if (nfgtasks.gt.1) then +c write (iout,*) "logPtot before reduction",logPtot + call MPI_Reduce(logPtot,logPtot_,1,MPI_DOUBLE_PRECISION, + & MPI_SUM,king,FG_COMM,IERR) + logPtot = logPtot_ +c write (iout,*) "logPtot after reduction",logPtot + call MPI_Reduce(gsaxsC(1,1),gsaxsC_(1,1),3*nres, + & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) + if (fg_rank.eq.king) then + do i=1,nres + do l=1,3 + gsaxsC(l,i) = gsaxsC_(l,i) + enddo + enddo + endif + call MPI_Reduce(gsaxsX(1,1),gsaxsX_(1,1),3*nres, + & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) + if (fg_rank.eq.king) then + do i=1,nres + do l=1,3 + gsaxsX(l,i) = gsaxsX_(l,i) + enddo + enddo + endif + endif +#endif + Esaxs_constr = logPtot + return + end