From: Adam Sieradzan Date: Tue, 29 Sep 2015 17:54:56 +0000 (+0200) Subject: Merge branch 'lipid' into AFM X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=25618f9f83673a7063414fe1e17415d138f58da8;hp=-c Merge branch 'lipid' into AFM Conflicts: source/wham/src-M/DIMENSIONS.ZSCOPT source/wham/src-M/energy_p_new.F source/wham/src-M/initialize_p.F source/wham/src-M/parmread.F --- 25618f9f83673a7063414fe1e17415d138f58da8 diff --combined source/unres/src_MD-M/energy_p_new_barrier.F index 5402971,3e1eb6e..ee55c93 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@@ -137,14 -137,6 +137,14 @@@ c print *,"Processor",myrank," com #ifdef TIMING time_vec=time_vec+MPI_Wtime()-time01 #endif +C Introduction of shielding effect first for each peptide group +C the shielding factor is set this factor is describing how each +C peptide group is shielded by side-chains +C the matrix - shield_fac(i) the i index describe the ith between i and i+1 +C write (iout,*) "shield_mode",shield_mode + if (shield_mode.gt.0) then + call set_shield_fac + endif c print *,"Processor",myrank," left VEC_AND_DERIV" if (ipot.lt.6) then #ifdef SPLITELE @@@ -201,10 -193,9 +201,10 @@@ C Calculate the virtual-bond-angle energy. C if (wang.gt.0d0) then - call ebend(ebe) + call ebend(ebe,ethetacnstr) else ebe=0 + ethetacnstr=0 endif c print *,"Processor",myrank," computed UB" C @@@ -330,7 -321,6 +330,7 @@@ energia(21)=esccor energia(22)=eliptran energia(23)=Eafmforce + energia(24)=ethetacnstr 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" @@@ -424,7 -414,6 +424,7 @@@ cMS$ATTRIBUTES C :: proc_pro esccor=energia(21) eliptran=energia(22) Eafmforce=energia(23) + ethetacnstr=energia(24) #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc @@@ -432,7 -421,6 +432,7 @@@ & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce + & +ethetacnstr #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc @@@ -441,7 -429,6 +441,7 @@@ & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran & +Eafmforce + & +ethetacnstr #endif energia(0)=etot c detecting NaNQ @@@ -546,7 -533,6 +546,7 @@@ c endd & wstrain*ghpbc(j,i) & +wliptran*gliptranc(j,i) & +gradafm(j,i) + & +welec*gshieldc(j,i) enddo enddo @@@ -565,7 -551,6 +565,7 @@@ & wstrain*ghpbc(j,i) & +wliptran*gliptranc(j,i) & +gradafm(j,i) + & +welec*gshieldc(j,i) enddo enddo @@@ -684,13 -669,6 +684,13 @@@ c endd do i=-1,nct do j=1,3 #ifdef SPLITELE +C print *,gradbufc(1,13) +C print *,welec*gelc(1,13) +C print *,wel_loc*gel_loc(1,13) +C print *,0.5d0*(wscp*gvdwc_scpp(1,13)) +C print *,welec*gelc_long(1,13)+wvdwpp*gvdwpp(1,13) +C print *,wel_loc*gel_loc_long(1,13) +C print *,gradafm(1,13),"AFM" gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & wel_loc*gel_loc(j,i)+ & 0.5d0*(wscp*gvdwc_scpp(j,i)+ @@@ -711,10 -689,6 +711,10 @@@ & +wscloc*gscloc(j,i) & +wliptran*gliptranc(j,i) & +gradafm(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + + #else gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & wel_loc*gel_loc(j,i)+ @@@ -736,9 -710,6 +736,9 @@@ & +wscloc*gscloc(j,i) & +wliptran*gliptranc(j,i) & +gradafm(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + #endif gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ @@@ -747,7 -718,6 +747,7 @@@ & wsccor*gsccorx(j,i) & +wscloc*gsclocx(j,i) & +wliptran*gliptranx(j,i) + & +welec*gshieldx(j,i) enddo enddo #ifdef DEBUG @@@ -1038,16 -1008,15 +1038,16 @@@ C-------------------------------------- esccor=energia(21) eliptran=energia(22) Eafmforce=energia(23) + ethetacnstr=energia(24) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp, & estr,wbond,ebe,wang, & escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain, & ecorr,wcorr, & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, - & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor, - & edihcnstr,ebr*nss, - & Uconst,eliptran,wliptran,Eafmforce,etot + & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr, + & ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, + & etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@@ -1069,7 -1038,6 +1069,7 @@@ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST= ',1pE16.6,' (Constraint energy)'/ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ @@@ -1083,8 -1051,7 +1083,8 @@@ & ecorr,wcorr, & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr, - & ebr*nss,Uconst,eliptran,wliptran,Eafmforc,etot + & ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, + & etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@@ -1105,7 -1072,6 +1105,7 @@@ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST=',1pE16.6,' (Constraint energy)'/ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ @@@ -1502,7 -1468,6 +1502,7 @@@ include 'COMMON.SBRIDGE' logical lprn integer xshift,yshift,zshift + evdw=0.0D0 ccccc energy_dec=.false. C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@@ -1599,36 -1564,10 +1599,36 @@@ do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN + +c write(iout,*) "PRZED ZWYKLE", evdwij call dyn_ssbond_ene(i,j,evdwij) +c write(iout,*) "PO ZWYKLE", evdwij + evdw=evdw+evdwij if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') & 'evdw',i,j,evdwij,' ss' +C triple bond artifac removal + do k=j+1,iend(i,iint) +C search over all next residues + if (dyn_ss_mask(k)) then +C check if they are cysteins +C write(iout,*) 'k=',k + +c write(iout,*) "PRZED TRI", evdwij + evdwij_przed_tri=evdwij + call triple_ssbond_ene(i,j,k,evdwij) +c if(evdwij_przed_tri.ne.evdwij) then +c write (iout,*) "TRI:", evdwij, evdwij_przed_tri +c endif + +c write(iout,*) "PO TRI", evdwij +C call the energy function that removes the artifical triple disulfide +C bond the soubroutine is located in ssMD.F + evdw=evdw+evdwij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') + & 'evdw',i,j,evdwij,'tss' + endif!dyn_ss_mask(k) + enddo! k ELSE ind=ind+1 itypj=iabs(itype(j)) @@@ -2780,6 -2719,17 +2780,17 @@@ c write(iout,*) 'b1=',b1(1,i-2 c write (iout,*) 'theta=', theta(i-1) enddo #else + if (i.gt. nnt+2 .and. i.lt.nct+2) then + iti = itortyp(itype(i-2)) + else + iti=ntortyp+1 + endif + c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then + if (i.gt. nnt+1 .and. i.lt.nct+1) then + iti1 = itortyp(itype(i-1)) + else + iti1=ntortyp+1 + endif b1(1,i-2)=b(3,iti) b1(2,i-2)=b(5,iti) b2(1,i-2)=b(2,iti) @@@ -2934,6 -2884,7 +2945,7 @@@ c if (i.gt. iatel_s+1 .and. i.lt do k=1,2 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,*) 'mu1',mu1(:,i-2) cd write (iout,*) 'mu2',mu2(:,i-2) @@@ -3446,9 -3397,7 +3458,9 @@@ C do zshift=-1, c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c +CTU KURWA do i=iatel_s,iatel_e +C do i=75,75 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 @@@ -3505,9 -3454,7 +3517,9 @@@ c endi c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) +C I TU KURWA do j=ielstart(i),ielend(i) +C do j=16,17 C write (iout,*) i,j if (j.le.1) cycle if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1 @@@ -3557,7 -3504,6 +3569,7 @@@ C-------------------------------------- include 'COMMON.FFIELD' include 'COMMON.TIME1' include 'COMMON.SPLITELE' + include 'COMMON.SHIELD' dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), @@@ -3690,22 -3636,10 +3702,22 @@@ c 4/26/02 - AL scaling down 1,4 repulsi el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)) el2=fac4*fac C MARYSIA - eesij=(el1+el2) +C eesij=(el1+el2) C 12/26/95 - for the evaluation of multi-body H-bonding interactions ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg) + if (shield_mode.gt.0) then +C fac_shield(i)=0.4 +C fac_shield(j)=0.6 + el1=el1*fac_shield(i)*fac_shield(j) + el2=el2*fac_shield(i)*fac_shield(j) + eesij=(el1+el2) + ees=ees+eesij + else + fac_shield(i)=1.0 + fac_shield(j)=1.0 + eesij=(el1+el2) ees=ees+eesij + endif evdw1=evdw1+evdwij*sss cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, @@@ -3729,99 -3663,22 +3741,99 @@@ erij(1)=xj*rmij erij(2)=yj*rmij erij(3)=zj*rmij + * * Radial derivatives. First process both termini of the fragment (i,j) * ggg(1)=facel*xj ggg(2)=facel*yj ggg(3)=facel*zj + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. + & (shield_mode.gt.0)) then +C print *,i,j + do ilist=1,ishield_list(i) + iresshield=shield_list(ilist,i) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i) + gshieldx(k,iresshield)=gshieldx(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) + gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield +C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield) +C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) +C if (iresshield.gt.i) then +C do ishi=i+1,iresshield-1 +C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield +C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) +C +C enddo +C else +C do ishi=iresshield,i +C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield +C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) +C +C enddo +C endif + enddo + enddo + do ilist=1,ishield_list(j) + iresshield=shield_list(ilist,j) + do k=1,3 + rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j) + gshieldx(k,iresshield)=gshieldx(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) + gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield + +C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) +C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield) +C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) +C if (iresshield.gt.j) then +C do ishi=j+1,iresshield-1 +C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield +C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) +C +C enddo +C else +C do ishi=iresshield,j +C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield +C & -grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) +C enddo +C endif + enddo + enddo + + do k=1,3 + gshieldc(k,i)=gshieldc(k,i)+ + & grad_shield(k,i)*eesij/fac_shield(i) + gshieldc(k,j)=gshieldc(k,j)+ + & grad_shield(k,j)*eesij/fac_shield(j) + gshieldc(k,i-1)=gshieldc(k,i-1)+ + & grad_shield(k,i)*eesij/fac_shield(i) + gshieldc(k,j-1)=gshieldc(k,j-1)+ + & grad_shield(k,j)*eesij/fac_shield(j) + + enddo + endif c do k=1,3 c ghalf=0.5D0*ggg(k) c gelc(k,i)=gelc(k,i)+ghalf c gelc(k,j)=gelc(k,j)+ghalf c enddo c 9/28/08 AL Gradient compotents will be summed only at the end +C print *,"before", gelc_long(1,i), gelc_long(1,j) do k=1,3 gelc_long(k,j)=gelc_long(k,j)+ggg(k) +C & +grad_shield(k,j)*eesij/fac_shield(j) gelc_long(k,i)=gelc_long(k,i)-ggg(k) +C & +grad_shield(k,i)*eesij/fac_shield(i) +C gelc_long(k,i-1)=gelc_long(k,i-1) +C & +grad_shield(k,i)*eesij/fac_shield(i) +C gelc_long(k,j-1)=gelc_long(k,j-1) +C & +grad_shield(k,j)*eesij/fac_shield(j) enddo +C print *,"bafter", gelc_long(1,i), gelc_long(1,j) + * * Loop over residues i+1 thru j-1. * @@@ -3870,11 -3727,8 +3882,11 @@@ C MARYSI * Radial derivatives. First process both termini of the fragment (i,j) * ggg(1)=fac*xj +C+eesij*grad_shield(1,i)+eesij*grad_shield(1,j) ggg(2)=fac*yj +C+eesij*grad_shield(2,i)+eesij*grad_shield(2,j) ggg(3)=fac*zj +C+eesij*grad_shield(3,i)+eesij*grad_shield(3,j) c do k=1,3 c ghalf=0.5D0*ggg(k) c gelc(k,i)=gelc(k,i)+ghalf @@@ -3917,8 -3771,7 +3929,8 @@@ c 9/28/08 AL Gradient compotents will b cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), cd & (dcosg(k),k=1,3) do k=1,3 - ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k) + ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))* + & fac_shield(i)*fac_shield(j) enddo c do k=1,3 c ghalf=0.5D0*ggg(k) @@@ -3934,21 -3787,16 +3946,21 @@@ cgrad do l=1, cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo +C print *,"before22", gelc_long(1,i), gelc_long(1,j) do k=1,3 gelc(k,i)=gelc(k,i) - & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) - & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) + & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) + & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)) + & *fac_shield(i)*fac_shield(j) gelc(k,j)=gelc(k,j) - & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) - & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) + & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) + & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)) + & *fac_shield(i)*fac_shield(j) gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo +C print *,"before33", gelc_long(1,i), gelc_long(1,j) + C MARYSIA c endif !sscale IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 @@@ -4151,7 -3999,7 +4163,7 @@@ C Contribution to the local-electrostat & +a33*muij(4) c write (iout,*) 'i',i,' j',j,itype(i),itype(j), c & ' eel_loc_ij',eel_loc_ij - c write(iout,*) 'muije=',muij(1),muij(2),muij(3),muij(4) + C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4) C Calculate patrial derivative for theta angle #ifdef NEWCORR geel_loc_ij=a22*gmuij1(1) @@@ -5273,13 -5121,8 +5285,13 @@@ include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' + include 'COMMON.CONTROL' dimension ggg(3) ehpb=0.0D0 + do i=1,3 + ggg(i)=0.0d0 + enddo +C write (iout,*) ,"link_end",link_end,constr_dist cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr cd write(iout,*)'link_start=',link_start,' link_end=',link_end if (link_end.eq.0) return @@@ -5306,84 -5149,27 +5318,84 @@@ cmc if (ii.gt.nres .and. itype(i C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds if (.not.dyn_ss .and. i.le.nss) then C 15/02/13 CC dynamic SSbond - additional check - if (ii.gt.nres - & .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then + if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. + & iabs(itype(jjj)).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij endif cd write (iout,*) "eij",eij +cd & ' waga=',waga,' fac=',fac + else if (ii.gt.nres .and. jj.gt.nres) then +c Restraints from contact prediction + dd=dist(ii,jj) + if (constr_dist.eq.11) then + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd + if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, + & ehpb,fordepth(i),dd + else + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "beta nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else + dd=dist(ii,jj) + rdis=dd-dhpb(i) +C Get the force constant corresponding to this distance. + waga=forcon(i) +C Calculate the contribution to energy. + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "beta reg",dd,waga*rdis*rdis +C +C Evaluate gradient. +C + fac=waga*rdis/dd + endif + endif + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo + do j=1,3 + ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) + ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) + enddo + do k=1,3 + ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) + enddo else C Calculate the distance between the two points and its difference from the C target distance. dd=dist(ii,jj) + if (constr_dist.eq.11) then + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd + if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, + & ehpb,fordepth(i),dd + else + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "alph nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else rdis=dd-dhpb(i) C Get the force constant corresponding to this distance. waga=forcon(i) C Calculate the contribution to energy. ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "alpha reg",dd,waga*rdis*rdis C C Evaluate gradient. C fac=waga*rdis/dd -cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, -cd & ' waga=',waga,' fac=',fac + endif + endif do j=1,3 ggg(j)=fac*(c(j,jj)-c(j,ii)) enddo @@@ -5407,7 -5193,7 +5419,7 @@@ cgrad endd enddo endif enddo - ehpb=0.5D0*ehpb + if (constr_dist.ne.11) ehpb=0.5D0*ehpb return end C-------------------------------------------------------------------------- @@@ -5598,7 -5384,7 +5610,7 @@@ end #ifdef CRYST_THETA C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@@ -5615,7 -5401,6 +5627,7 @@@ include 'COMMON.NAMES' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.TORCNSTR' common /calcthet/ term1,term2,termm,diffak,ratak, & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq, & delthe0,sig0inv,sigtc,sigsqtc,delthec,it @@@ -5727,34 -5512,6 +5739,34 @@@ C Derivatives of the "mean" values in g if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg) enddo + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=ithetaconstr_start,ithetaconstr_end + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif + if (energy_dec) then + write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", + & i,itheta,rad2deg*thetiii, + & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), + & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, + & gloc(itheta+nphi-2,icg) + endif + enddo + C Ufff.... We've done all this!!! return end @@@ -5871,7 -5628,7 +5883,7 @@@ C "Thank you" to MAPLE (probably spare end #else C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@@ -5890,7 -5647,6 +5902,7 @@@ include 'COMMON.NAMES' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.TORCNSTR' double precision coskt(mmaxtheterm),sinkt(mmaxtheterm), & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle), & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble), @@@ -5901,7 -5657,8 +5913,7 @@@ 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 -C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF - +C print *,i,theta(i) if (iabs(itype(i+1)).eq.20) iblock=2 if (iabs(itype(i+1)).ne.20) iblock=1 dethetai=0.0d0 @@@ -5913,7 -5670,6 +5925,7 @@@ coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo +C print *,ethetai if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) @@@ -5929,8 -5685,8 +5941,8 @@@ C propagation of chirality for glycine enddo else phii=0.0d0 - ityp1=nthetyp+1 do k=1,nsingle + ityp1=ithetyp((itype(i-2))) cosph1(k)=0.0d0 sinph1(k)=0.0d0 enddo @@@ -5950,7 -5706,7 +5962,7 @@@ enddo else phii1=0.0d0 - ityp3=nthetyp+1 + ityp3=ithetyp((itype(i))) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@@ -6000,7 -5756,6 +6012,7 @@@ enddo write(iout,*) "ethetai",ethetai endif +C print *,ethetai do m=1,ntheterm2 do k=1,nsingle aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k) @@@ -6021,16 -5776,10 +6033,16 @@@ & ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet", & ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet", & eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai +C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k) enddo enddo +C print *,"cosph1", (cosph1(k), k=1,nsingle) +C print *,"cosph2", (cosph2(k), k=1,nsingle) +C print *,"sinph1", (sinph1(k), k=1,nsingle) +C print *,"sinph2", (sinph2(k), k=1,nsingle) if (lprn) & write(iout,*) "ethetai",ethetai +C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k) do m=1,ntheterm3 do k=2,ndouble do l=1,k-1 @@@ -6066,7 -5815,6 +6078,7 @@@ enddo 10 continue c lprn1=.true. +C print *,ethetai if (lprn1) & write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') & i,theta(i)*rad2deg,phii*rad2deg, @@@ -6075,37 -5823,8 +6087,37 @@@ c lprn1=.false etheta=etheta+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)=wang*dethetai+gloc(nphi+i-2,icg) + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai + enddo +C now constrains + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=ithetaconstr_start,ithetaconstr_end + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif + if (energy_dec) then + write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", + & i,itheta,rad2deg*thetiii, + & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), + & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, + & gloc(itheta+nphi-2,icg) + endif enddo + return end #endif @@@ -6925,12 -6644,12 +6937,12 @@@ c write (iout,*) 'i=',i,' gloc=', difi=phii-phi0(i) if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)**difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 endif ! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, ! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) @@@ -7036,21 -6755,18 +7048,21 @@@ c do i=1,ndih_const difi=pinorm(phii-phi0(i)) if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else difi=0.0 endif -cd write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii, -cd & rad2deg*phi0(i), rad2deg*drange(i), -cd & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) + if (energy_dec) then + write (iout,'(a6,2i5,4f8.3,2e14.5)') "edihc", + & i,itori,rad2deg*phii, + & rad2deg*phi0(i), rad2deg*drange(i), + & rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg) + endif enddo cd write (iout,*) 'edihcnstr',edihcnstr return @@@ -9127,9 -8843,9 +9139,9 @@@ cd ghalf=0.0d cold ghalf=0.5d0*eel5*eij*gacont_hbr(ll,kk,k) cgrad ghalf=0.5d0*ggg2(ll) cd ghalf=0.0d0 - gradcorr5(ll,k)=gradcorr5(ll,k)+ghalf+ekont*derx(ll,2,2) + gradcorr5(ll,k)=gradcorr5(ll,k)+ekont*derx(ll,2,2) gradcorr5(ll,k+1)=gradcorr5(ll,k+1)+ekont*derx(ll,3,2) - gradcorr5(ll,l)=gradcorr5(ll,l)+ghalf+ekont*derx(ll,4,2) + gradcorr5(ll,l)=gradcorr5(ll,l)+ekont*derx(ll,4,2) gradcorr5(ll,l1)=gradcorr5(ll,l1)+ekont*derx(ll,5,2) gradcorr5_long(ll,l)=gradcorr5_long(ll,l)+gradcorr5kl gradcorr5_long(ll,k)=gradcorr5_long(ll,k)-gradcorr5kl @@@ -10633,164 -10349,4 +10645,164 @@@ C Eafmforce=-forceAFMconst*(dist-d C print *,'AFM',Eafmforce,totTafm*velAFMconst,dist return end +C----------------------------------------------------------- +C first for shielding is setting of function of side-chains + subroutine set_shield_fac + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.IOUNITS' + include 'COMMON.SHIELD' + include 'COMMON.INTERACT' +C this is the squar root 77 devided by 81 the epislion in lipid (in protein) + double precision div77_81/0.974996043d0/, + &div4_81/0.2222222222d0/,sh_frac_dist_grad(3) + +C the vector between center of side_chain and peptide group + double precision pep_side(3),long,side_calf(3), + &pept_group(3),costhet_grad(3),cosphi_grad_long(3), + &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3) +C the line belowe needs to be changed for FGPROC>1 + do i=1,nres-1 + if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle + ishield_list(i)=0 +Cif there two consequtive dummy atoms there is no peptide group between them +C the line below has to be changed for FGPROC>1 + VolumeTotal=0.0 + do k=1,nres + if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle + dist_pep_side=0.0 + dist_side_calf=0.0 + do j=1,3 +C first lets set vector conecting the ithe side-chain with kth side-chain + pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0 +C pep_side(j)=2.0d0 +C and vector conecting the side-chain with its proper calfa + side_calf(j)=c(j,k+nres)-c(j,k) +C side_calf(j)=2.0d0 + pept_group(j)=c(j,i)-c(j,i+1) +C lets have their lenght + dist_pep_side=pep_side(j)**2+dist_pep_side + dist_side_calf=dist_side_calf+side_calf(j)**2 + dist_pept_group=dist_pept_group+pept_group(j)**2 + enddo + dist_pep_side=dsqrt(dist_pep_side) + dist_pept_group=dsqrt(dist_pept_group) + dist_side_calf=dsqrt(dist_side_calf) + do j=1,3 + pep_side_norm(j)=pep_side(j)/dist_pep_side + side_calf_norm(j)=dist_side_calf + enddo +C now sscale fraction + sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield +C print *,buff_shield,"buff" +C now sscale + if (sh_frac_dist.le.0.0) cycle +C If we reach here it means that this side chain reaches the shielding sphere +C Lets add him to the list for gradient + ishield_list(i)=ishield_list(i)+1 +C ishield_list is a list of non 0 side-chain that contribute to factor gradient +C this list is essential otherwise problem would be O3 + shield_list(ishield_list(i),i)=k +C Lets have the sscale value + if (sh_frac_dist.gt.1.0) then + scale_fac_dist=1.0d0 + do j=1,3 + sh_frac_dist_grad(j)=0.0d0 + enddo + else + scale_fac_dist=-sh_frac_dist*sh_frac_dist + & *(2.0*sh_frac_dist-3.0d0) + fac_help_scale=6.0*(sh_frac_dist-sh_frac_dist**2) + & /dist_pep_side/buff_shield*0.5 +C remember for the final gradient multiply sh_frac_dist_grad(j) +C for side_chain by factor -2 ! + do j=1,3 + sh_frac_dist_grad(j)=fac_help_scale*pep_side(j) +C print *,"jestem",scale_fac_dist,fac_help_scale, +C & sh_frac_dist_grad(j) + enddo + endif +C if ((i.eq.3).and.(k.eq.2)) then +C print *,i,sh_frac_dist,dist_pep,fac_help_scale,scale_fac_dist +C & ,"TU" +C endif + +C this is what is now we have the distance scaling now volume... + short=short_r_sidechain(itype(k)) + long=long_r_sidechain(itype(k)) + costhet=1.0d0/dsqrt(1.0+short**2/dist_pep_side**2) +C now costhet_grad +C costhet=0.0d0 + costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side**4 +C costhet_fac=0.0d0 + do j=1,3 + costhet_grad(j)=costhet_fac*pep_side(j) + enddo +C remember for the final gradient multiply costhet_grad(j) +C for side_chain by factor -2 ! +C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1 +C pep_side0pept_group is vector multiplication + pep_side0pept_group=0.0 + do j=1,3 + pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j) + enddo + cosalfa=(pep_side0pept_group/ + & (dist_pep_side*dist_side_calf)) + fac_alfa_sin=1.0-cosalfa**2 + fac_alfa_sin=dsqrt(fac_alfa_sin) + rkprim=fac_alfa_sin*(long-short)+short +C now costhet_grad + cosphi=1.0d0/dsqrt(1.0+rkprim**2/dist_pep_side**2) + cosphi_fac=cosphi**3*rkprim**2*(-0.5)/dist_pep_side**4 + + do j=1,3 + cosphi_grad_long(j)=cosphi_fac*pep_side(j) + &+cosphi**3*0.5/dist_pep_side**2*(-rkprim) + &*(long-short)/fac_alfa_sin*cosalfa/ + &((dist_pep_side*dist_side_calf))* + &((side_calf(j))-cosalfa* + &((pep_side(j)/dist_pep_side)*dist_side_calf)) + + cosphi_grad_loc(j)=cosphi**3*0.5/dist_pep_side**2*(-rkprim) + &*(long-short)/fac_alfa_sin*cosalfa + &/((dist_pep_side*dist_side_calf))* + &(pep_side(j)- + &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side) + enddo + + VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi) + & /VSolvSphere_div +C now the gradient... +C grad_shield is gradient of Calfa for peptide groups + do j=1,3 + grad_shield(j,i)=grad_shield(j,i) +C gradient po skalowaniu + & +(sh_frac_dist_grad(j) +C gradient po costhet + &-scale_fac_dist*costhet_grad(j)/(1.0-costhet) + &-scale_fac_dist*(cosphi_grad_long(j)) + &/(1.0-cosphi) )*div77_81 + &*VofOverlap +C grad_shield_side is Cbeta sidechain gradient + grad_shield_side(j,ishield_list(i),i)= + & (sh_frac_dist_grad(j)*-2.0d0 + & +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet) + & +scale_fac_dist*(cosphi_grad_long(j)) + & *2.0d0/(1.0-cosphi)) + & *div77_81*VofOverlap + + grad_shield_loc(j,ishield_list(i),i)= + & scale_fac_dist*cosphi_grad_loc(j) + & *2.0d0/(1.0-cosphi) + & *div77_81*VofOverlap + enddo + VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist + enddo + fac_shield(i)=VolumeTotal*div77_81+div4_81 +C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i) + enddo + return + end diff --combined source/unres/src_MD-M/parmread.F index c881425,183a917..ae4d710 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@@ -26,8 -26,6 +26,8 @@@ include 'COMMON.SBRIDGE' include 'COMMON.MD' include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.SHIELD' character*1 t1,t2,t3 character*1 onelett(4) /"G","A","P","D"/ character*1 toronelet(-2:2) /"p","a","G","A","P"/ @@@ -75,7 -73,6 +75,7 @@@ #else read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok do i=1,ntyp + print *,i read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i), & j=1,nbondterm(i)),msc(i),isc(i),restok(i) dsc(i) = vbldsc0(1,i) @@@ -947,16 -944,6 +947,16 @@@ c B2(1,i) = b(2 c B2(2,i) = b(4) c B2(1,-i) =b(2) c B2(2,-i) =-b(4) + B1tilde(1,i) = b(3,i) + B1tilde(2,i) =-b(5,i) +C B1tilde(1,-i) =-b(3,i) +C B1tilde(2,-i) =b(5,i) + b1tilde(1,i)=0.0d0 + b1tilde(2,i)=0.0d0 + B2(1,i) = b(2,i) + B2(2,i) = b(4,i) +C B2(1,-i) =b(2,i) +C B2(2,-i) =-b(4,i) c b2(1,i)=0.0d0 c b2(2,i)=0.0d0 @@@ -1095,7 -1082,7 +1095,7 @@@ & ', exponents are ',expon,2*expon goto (10,20,30,30,40) ipot C----------------------- LJ potential --------------------------------- - 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp), + 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp), & (sigma0(i),i=1,ntyp) if (lprint) then write (iout,'(/a/)') 'Parameters of the LJ potential:' @@@ -1107,7 -1094,7 +1107,7 @@@ endif goto 50 C----------------------- LJK potential -------------------------------- - 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp), + 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp), & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp) if (lprint) then write (iout,'(/a/)') 'Parameters of the LJK potential:' @@@ -1121,7 -1108,7 +1121,7 @@@ goto 50 C---------------------- GB or BP potential ----------------------------- 30 do i=1,ntyp - read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp) + read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp) enddo read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp) read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp) @@@ -1130,11 -1117,12 +1130,12 @@@ C now we start reading lipid do i=1,ntyp read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp) - print *,"WARNING!!" - do j=1,ntyp - epslip(i,j)=epslip(i,j)+0.05d0 - enddo + C print *,"WARNING!!" + C do j=1,ntyp + C epslip(i,j)=epslip(i,j)+0.05d0 + C enddo enddo + write(iout,*) epslip(1,1),"OK?" C For the GB potential convert sigma'**2 into chi' if (ipot.eq.4) then do i=1,ntyp @@@ -1153,7 -1141,7 +1154,7 @@@ endif goto 50 C--------------------- GBV potential ----------------------------------- - 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp), + 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp), & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp), & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp) if (lprint) then @@@ -1302,7 -1290,7 +1303,7 @@@ c lprint=.false C C Define the constants of the disulfide bridge C - ebr=-5.50D0 +C ebr=-12.00D0 c c Old arbitrary potential - commented out. c @@@ -1313,13 -1301,13 +1314,13 @@@ c Constants of the disulfide-bond poten c energy surface of diethyl disulfide. c A. Liwo and U. Kozlowska, 11/24/03 c - D0CM = 3.78d0 - AKCM = 15.1d0 - AKTH = 11.0d0 - AKCT = 12.0d0 - V1SS =-1.08d0 - V2SS = 7.61d0 - V3SS = 13.7d0 +C D0CM = 3.78d0 +C AKCM = 15.1d0 +C AKTH = 11.0d0 +C AKCT = 12.0d0 +C V1SS =-1.08d0 +C V2SS = 7.61d0 +C V3SS = 13.7d0 c akcm=0.0d0 c akth=0.0d0 c akct=0.0d0 @@@ -1327,33 -1315,14 +1328,33 @@@ c v1ss=0.0d c v2ss=0.0d0 c v3ss=0.0d0 - if(me.eq.king) then - write (iout,'(/a)') "Disulfide bridge parameters:" - write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr - write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm - write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct - write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, - & ' v3ss:',v3ss - endif +C if(me.eq.king) then +C write (iout,'(/a)') "Disulfide bridge parameters:" +C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr +C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm +C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct +C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, +C & ' v3ss:',v3ss +C endif +C set the variables used for shielding effect +C write (iout,*) "SHIELD MODE",shield_mode +C if (shield_mode.gt.0) then +C VSolvSphere the volume of solving sphere +C print *,pi,"pi" +C rpp(1,1) is the energy r0 for peptide group contact and will be used for it +C there will be no distinction between proline peptide group and normal peptide +C group in case of shielding parameters +C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 +C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 +C write (iout,*) VSolvSphere,VSolvSphere_div +C long axis of side chain +C do i=1,ntyp +C long_r_sidechain(i)=vbldsc0(1,i) +C short_r_sidechain(i)=sigma0(i) +C enddo +C lets set the buffor value +C buff_shield=1.0d0 +C endif return 111 write (iout,*) "Error reading bending energy parameters." goto 999 @@@ -1425,22 -1394,6 +1426,22 @@@ c-HP- if(ierror.ne.0) stop '--error ret #else call getenv(var,val) #endif - +C set the variables used for shielding effect +C if (shield_mode.gt.0) then +C VSolvSphere the volume of solving sphere +C print *,pi,"pi" +C rpp(1,1) is the energy r0 for peptide group contact and will be used for it +C there will be no distinction between proline peptide group and normal peptide +C group in case of shielding parameters +C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 +C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 +C long axis of side chain +C do i=1,ntyp +C long_r_sidechain(i)=vbldsc0(1,i) +C short_r_sidechain(i)=sigma0(i) +C enddo +C lets set the buffor value +C buff_shield=1.0d0 +C endif return end diff --combined source/unres/src_MD-M/readrtns_CSA.F index d25c4eb,889e5a5..ebb59da --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@@ -81,7 -81,6 +81,7 @@@ include 'COMMON.INTERACT' include 'COMMON.SETUP' include 'COMMON.SPLITELE' + include 'COMMON.SHIELD' COMMON /MACHSW/ KDIAG,ICORFL,IXDR character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/ character*80 ucase @@@ -99,10 -98,6 +99,10 @@@ c print *,"Processor",me," fg_rank C Set up the time limit (caution! The time must be input in minutes!) read_cart=index(controlcard,'READ_CART').gt.0 call readi(controlcard,'CONSTR_DIST',constr_dist,0) +C this variable with_theta_constr is the variable which allow to read and execute the +C constrains on theta angles WITH_THETA_CONSTR is the keyword + with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0 + write (iout,*) "with_theta_constr ",with_theta_constr call readi(controlcard,'SYM',symetr,1) call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0 @@@ -146,15 -141,6 +146,15 @@@ selfguide=(index(controlcard,'SELFGUIDE')) print *,'AFMlog',AFMlog,selfguide,"KUPA" call readi(controlcard,'IPRINT',iprint,0) +C SHIELD keyword sets if the shielding effect of side-chains is used +C 0 denots no shielding is used all peptide are equally despite the +C solvent accesible area +C 1 the newly introduced function +C 2 reseved for further possible developement + call readi(controlcard,'SHIELD',shield_mode,0) +C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then + write(iout,*) "shield_mode",shield_mode +C endif call readi(controlcard,'MAXGEN',maxgen,10000) call readi(controlcard,'MAXOVERLAP',maxoverlap,1000) call readi(controlcard,"KDIAG",kdiag,0) @@@ -247,7 -233,7 +247,7 @@@ c Cutoff range for interaction bordliptop=(boxzsize+lipthick)/2.0 bordlipbot=bordliptop-lipthick C endif - if ((bordliptop.gt.boxzsize).or.(borlipbot.lt.0.0)) + if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" buflipbot=bordlipbot+lipbufthick bufliptop=bordliptop-lipbufthick @@@ -258,24 -244,8 +258,24 @@@ write(iout,*) "bordlipbot=",bordlipbot write(iout,*) "bufliptop=",bufliptop write(iout,*) "buflipbot=",buflipbot - - + write (iout,*) "SHIELD MODE",shield_mode + if (shield_mode.gt.0) then + pi=3.141592d0 +C VSolvSphere the volume of solving sphere +C print *,pi,"pi" +C rpp(1,1) is the energy r0 for peptide group contact and will be used for it +C there will be no distinction between proline peptide group and normal peptide +C group in case of shielding parameters + VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 + VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 + write (iout,*) VSolvSphere,VSolvSphere_div +C long axis of side chain + do i=1,ntyp + long_r_sidechain(i)=vbldsc0(1,i) + short_r_sidechain(i)=sigma0(i) + enddo + buff_shield=1.0d0 + endif if (me.eq.king .or. .not.out1file ) & write (iout,*) "DISTCHAINMAX",distchainmax @@@ -601,7 -571,7 +601,7 @@@ integer rescode double precision x(maxvar) character*256 pdbfile - character*320 weightcard + character*400 weightcard character*80 weightcard_t,ucase dimension itype_pdb(maxres) common /pizda/ itype_pdb @@@ -743,14 -713,6 +743,14 @@@ C 12/1/95 Added weight for the multi-bo call reada(weightcard,"V2SS",v2ss,7.61d0) call reada(weightcard,"V3SS",v3ss,13.7d0) call reada(weightcard,"EBR",ebr,-5.50D0) + call reada(weightcard,"ATRISS",atriss,0.301D0) + call reada(weightcard,"BTRISS",btriss,0.021D0) + call reada(weightcard,"CTRISS",ctriss,1.001D0) + call reada(weightcard,"DTRISS",dtriss,1.001D0) + write (iout,*) "ATRISS=", atriss + write (iout,*) "BTRISS=", btriss + write (iout,*) "CTRISS=", ctriss + write (iout,*) "DTRISS=", dtriss dyn_ss=(index(weightcard,'DYN_SS').gt.0) do i=1,maxres dyn_ss_mask(i)=.false. @@@ -771,11 -733,7 +771,11 @@@ v2ss=v2ss*wstrain/wsc v3ss=v3ss*wstrain/wsc else - ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain + if (wstrain.ne.0.0) then + ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain + else + ss_depth=0.0 + endif endif if(me.eq.king.or..not.out1file) then @@@ -797,11 -755,9 +797,11 @@@ 33 write (iout,'(a)') 'Error opening PDB file.' stop 34 continue -c print *,'Begin reading pdb data' +c write (iout,*) 'Begin reading pdb data' +c call flush(iout) call readpdb -c print *,'Finished reading pdb data' +c write (iout,*) 'Finished reading pdb data' +c call flush(iout) if(me.eq.king.or..not.out1file) & write (iout,'(a,i3,a,i3)')'nsup=',nsup, & ' nstart_sup=',nstart_sup @@@ -891,78 -847,27 +891,78 @@@ C 8/13/98 Set limits to generating the enddo read (inp,*) ndih_constr if (ndih_constr.gt.0) then - read (inp,*) ftors - read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr) +C read (inp,*) ftors + read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i), + & i=1,ndih_constr) if(me.eq.king.or..not.out1file)then write (iout,*) & 'There are',ndih_constr,' constraints on phi angles.' do i=1,ndih_constr - write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i) + write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i), + & ftors(i) enddo endif do i=1,ndih_constr phi0(i)=deg2rad*phi0(i) drange(i)=deg2rad*drange(i) enddo - if(me.eq.king.or..not.out1file) - & write (iout,*) 'FTORS',ftors +C if(me.eq.king.or..not.out1file) +C & write (iout,*) 'FTORS',ftors do i=1,ndih_constr ii = idih_constr(i) phibound(1,ii) = phi0(i)-drange(i) phibound(2,ii) = phi0(i)+drange(i) enddo endif +C first setting the theta boundaries to 0 to pi +C this mean that there is no energy penalty for any angle occuring this can be applied +C for generate random conformation but is not implemented in this way +C do i=1,nres +C thetabound(1,i)=0 +C thetabound(2,i)=pi +C enddo +C begin reading theta constrains this is quartic constrains allowing to +C have smooth second derivative + if (with_theta_constr) then +C with_theta_constr is keyword allowing for occurance of theta constrains + read (inp,*) ntheta_constr +C ntheta_constr is the number of theta constrains + if (ntheta_constr.gt.0) then +C read (inp,*) ftors + read (inp,*) (itheta_constr(i),theta_constr0(i), + & theta_drange(i),for_thet_constr(i), + & i=1,ntheta_constr) +C the above code reads from 1 to ntheta_constr +C itheta_constr(i) residue i for which is theta_constr +C theta_constr0 the global minimum value +C theta_drange is range for which there is no energy penalty +C for_thet_constr is the force constant for quartic energy penalty +C E=k*x**4 + if(me.eq.king.or..not.out1file)then + write (iout,*) + & 'There are',ntheta_constr,' constraints on phi angles.' + do i=1,ntheta_constr + write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), + & theta_drange(i), + & for_thet_constr(i) + enddo + endif + do i=1,ntheta_constr + theta_constr0(i)=deg2rad*theta_constr0(i) + theta_drange(i)=deg2rad*theta_drange(i) + enddo +C if(me.eq.king.or..not.out1file) +C & write (iout,*) 'FTORS',ftors +C do i=1,ntheta_constr +C ii = itheta_constr(i) +C thetabound(1,ii) = phi0(i)-drange(i) +C thetabound(2,ii) = phi0(i)+drange(i) +C enddo + endif ! ntheta_constr.gt.0 + endif! with_theta_constr +C +C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0 +C write (iout,*) "with_dihed_constr ",with_dihed_constr nnt=1 #ifdef MPI if (me.eq.king) then @@@ -1049,9 -954,7 +1049,9 @@@ czscore call geom_to_var(nvar, enddo call contact(.true.,ncont_ref,icont_ref,co) endif -c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup + endif + print *, "A TU" + write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup call flush(iout) if (constr_dist.gt.0) call read_dist_constr write (iout,*) "After read_dist_constr nhpb",nhpb @@@ -1072,7 -975,7 +1072,7 @@@ & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i) enddo endif - endif +C endif if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and. & modecalc.ne.10) then @@@ -1132,7 -1035,6 +1132,7 @@@ C initial geometry omeg(i)=-120d0*deg2rad if (itype(i).le.0) omeg(i)=-omeg(i) enddo + call chainbuild_extconf else if(me.eq.king.or..not.out1file) & write (iout,'(a)') 'Random-generated initial geometry.' @@@ -1245,7 -1147,6 +1245,7 @@@ cd write (iout,'(i4,f10.5)') (i,rad2 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') & 'Processor',myrank,': end reading molecular data.' #endif + print *,"A TU?" return end c-------------------------------------------------------------------------- @@@ -2402,8 -2303,7 +2402,8 @@@ c-------------------------------------- integer ifrag_(2,100),ipair_(2,100) double precision wfrag_(100),wpair_(100) character*500 controlcard -c write (iout,*) "Calling read_dist_constr" + print *, "WCHODZE" + write (iout,*) "Calling read_dist_constr" c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup c call flush(iout) call card_concat(controlcard) @@@ -2497,30 -2397,11 +2497,30 @@@ c write (iout,*) "j",j," k", enddo endif enddo + print *,ndist_ do i=1,ndist_ - read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) + if (constr_dist.eq.11) then + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), + & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1) + fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) + else +C print *,"in else" + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), + & ibecarb(i),forcon(nhpb+1) + endif if (forcon(nhpb+1).gt.0.0d0) then nhpb=nhpb+1 - dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + if (ibecarb(i).gt.0) then + ihpb(i)=ihpb(i)+nres + jhpb(i)=jhpb(i)+nres + endif + if (dhpb(nhpb).eq.0.0d0) + & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + endif +C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) +C if (forcon(nhpb+1).gt.0.0d0) then +C nhpb=nhpb+1 +C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", @@@ -2529,7 -2410,7 +2529,7 @@@ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif - endif + enddo call flush(iout) return diff --combined source/wham/src-M/enecalc1.F index fee94cf,81668b5..a5d25b3 --- a/source/wham/src-M/enecalc1.F +++ b/source/wham/src-M/enecalc1.F @@@ -35,7 -35,7 +35,7 @@@ double precision tole /1.0d-1/ integer i,itj,ii,iii,j,k,l,licz integer ir,ib,ipar,iparm - integer iscor,islice + integer iscor,islice,scount_buff(0:99) real*4 csingle(3,maxres2) double precision energ double precision temp @@@ -167,7 -167,7 +167,7 @@@ C write (iout,*) "tuz za energia write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) call enerprint(energia(0),fT) write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,21) - write (iout,*) "ftors",ftors + write (iout,*) "ftors(1)",ftors(1) call briefout(i,energia(0)) temp=1.0d0/(beta_h(ib,ipar)*1.987D-3) write (iout,*) "temp", temp @@@ -228,7 -228,7 +228,7 @@@ c call intou endif C write (iout,*) "Czy tu dochodze" potE(iii+1,iparm)=energia(0) - do k=1,21 + do k=1,22 enetb(k,iii+1,iparm)=energia(k) enddo #ifdef DEBUG @@@ -277,15 -277,12 +277,15 @@@ c & " snk",snk_p(iR,ib,ipar 121 continue enddo #ifdef MPI - scount(me)=iii - write (iout,*) "Me",me," scount",scount(me) + scount_buff(me)=iii + write (iout,*) "Me",me," scount_buff",scount_buff(me) call flush(iout) c Master gathers updated numbers of conformations written by all procs. - call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount(0), 1, +c call MPI_AllGather(MPI_IN_PLACE,1,MPI_DATATYPE_NULL,scount(0),1, +c & MPI_INTEGER, WHAM_COMM, IERROR) + call MPI_AllGather( scount_buff(me), 1, MPI_INTEGER, scount(0), 1, & MPI_INTEGER, WHAM_COMM, IERROR) + indstart(0)=1 indend(0)=scount(0) do i=1, Nprocs-1 @@@ -371,7 -368,7 +371,7 @@@ c-------------------------------------- double precision energ integer ilen,iroof external ilen,iroof - integer ir,ib,iparm + integer ir,ib,iparm, scount_buff(0:99) integer isecstr(maxres) write (licz2,'(bz,i2.2)') islice call opentmp(islice,ientout,bprotfile_temp) @@@ -676,13 -673,8 +676,13 @@@ c write (iout,*) "xdrf3dfcoord c call flush(iout) call xdrfint_(ixdrf, nss, iret) do j=1,nss - call xdrfint_(ixdrf, ihpb(j), iret) - call xdrfint_(ixdrf, jhpb(j), iret) + if (dyn_ss) then + call xdrfint(ixdrf, idssb(j)+nres, iret) + call xdrfint(ixdrf, jdssb(j)+nres, iret) + else + call xdrfint_(ixdrf, ihpb(j), iret) + call xdrfint_(ixdrf, jhpb(j), iret) + endif enddo call xdrffloat_(ixdrf,real(eini),iret) call xdrffloat_(ixdrf,real(efree),iret) @@@ -693,13 -685,8 +693,13 @@@ call xdrfint(ixdrf, nss, iret) do j=1,nss - call xdrfint(ixdrf, ihpb(j), iret) - call xdrfint(ixdrf, jhpb(j), iret) + if (dyn_ss) then + call xdrfint(ixdrf, idssb(j)+nres, iret) + call xdrfint(ixdrf, jdssb(j)+nres, iret) + else + call xdrfint(ixdrf, ihpb(j), iret) + call xdrfint(ixdrf, jhpb(j), iret) + endif enddo call xdrffloat(ixdrf,real(eini),iret) call xdrffloat(ixdrf,real(efree),iret) diff --combined source/wham/src-M/energy_p_new.F index a7b0798,b8e07c5..cba6b5e --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@@ -70,9 -70,8 +70,9 @@@ cd print *,'EHPB exitted succesfully C C Calculate the virtual-bond-angle energy. C - call ebend(ebe) C print *,'Bend energy finished.' + call ebend(ebe,ethetacnstr) +cd print *,'Bend energy finished.' C C Calculate the SC local energy. C @@@ -91,6 -90,11 +91,11 @@@ C 21/5/07 Calculate local sicdechain correlation energy C call eback_sc_corr(esccor) + + if (wliptran.gt.0) then + call Eliptransfer(eliptran) + endif + C C 12/1/95 Multi-body terms C @@@ -111,20 -115,20 +116,22 @@@ c write (iout,*) "ft(6)",fact(6), etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees & +wvdwpp*evdw1 & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 + & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d - & +wbond*estr+wsccor*fact(1)*esccor+wliptran*eliptran + & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr ++ & +wliptran*eliptran #else etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2 & +welec*fact(1)*(ees+evdw1) & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 + & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d - & +wbond*estr+wsccor*fact(1)*esccor+wliptran*eliptran + & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr ++ & +wliptran*eliptran #endif energia(0)=etot energia(1)=evdw @@@ -158,7 -162,8 +165,8 @@@ energia(19)=esccor energia(20)=edihcnstr energia(21)=evdw_t + energia(24)=ethetacnstr + energia(22)=eliptran - c detecting NaNQ #ifdef ISNAN #ifdef AIX @@@ -197,10 -202,12 +205,12 @@@ & wcorr6*fact(5)*gradcorr6(j,i)+ & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) 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*fact(2)*gsccorx(j,i) + & +wliptran*gliptranx(j,i) enddo #else do i=1,nct @@@ -216,10 -223,12 +226,12 @@@ & wcorr6*fact(5)*gradcorr6(j,i)+ & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) + & +wliptran*gliptranc(j,i) 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*fact(1)*gsccorx(j,i) + & +wliptran*gliptranx(j,i) enddo #endif enddo @@@ -233,11 -242,8 +245,11 @@@ & +wturn3*fact(2)*gel_loc_turn3(i) & +wturn6*fact(5)*gel_loc_turn6(i) & +wel_loc*fact(2)*gel_loc_loc(i) +c & +wsccor*fact(1)*gsccor_loc(i) +c BYLA ROZNICA Z CLUSTER< OSTATNIA LINIA DODANA enddo endif + if (dyn_ss) call dyn_set_nss return end C------------------------------------------------------------------------ @@@ -275,7 -281,7 +287,8 @@@ esccor=energia(19) edihcnstr=energia(20) estr=energia(18) + ethetacnstr=energia(24) + eliptran=energia(22) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1, & wvdwpp, @@@ -284,7 -290,7 +297,8 @@@ & ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5), & eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2), & eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5), - & esccor,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss,etot - & esccor,wsccor*fact(1),edihcnstr,ebr*nss,eliptran,wliptran,etot ++ & esccor,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss, ++ & eliptran,wliptran,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@@ -306,8 -312,8 +320,9 @@@ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ + & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ & 'ETOT= ',1pE16.6,' (total)') #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),estr,wbond, @@@ -316,7 -322,7 +331,7 @@@ & ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2), & eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3), & eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor, - & edihcnstr,ethetacnstr,ebr*nss,etot - & edihcnstr,ebr*nss,eliptran,wliptran,etot ++ & edihcnstr,ethetacnstr,ebr*nss,eliptran,wliptran,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@@ -337,8 -343,8 +352,9 @@@ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ + & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@@ -370,14 -376,11 +386,14 @@@ integer icant external icant cd print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon +c ROZNICA z cluster do i=1,210 do j=1,2 eneps_temp(j,i)=0.0d0 enddo enddo +cROZNICA + evdw=0.0D0 evdw_t=0.0d0 do i=iatsc_s,iatsc_e @@@ -407,22 -410,19 +423,22 @@@ C Change 12/1/95 to calculate four-bod c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj eps0ij=eps(itypi,itypj) fac=rrij**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=e1+e2 ij=icant(itypi,itypj) +c ROZNICA z cluster eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij) eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps0ij +c + cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj) cd write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)') cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm, cd & (c(k,i),k=1,3),(c(k,j),k=1,3) - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij else evdw_t=evdw_t+evdwij @@@ -580,8 -580,8 +596,8 @@@ rij=1.0D0/r_inv_ij r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=e_augm+e1+e2 ij=icant(itypi,itypj) eneps_temp(1,ij)=eneps_temp(1,ij)+(e1+a_augm) @@@ -594,7 -594,7 +610,7 @@@ cd & restyp(itypi),i,restyp(it cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm, cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij, cd & (c(k,i),k=1,3),(c(k,j),k=1,3) - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij else evdw_t=evdw_t+evdwij @@@ -726,8 -726,8 +742,8 @@@ C Calculate the angle-dependent terms o C Calculate whole angle-dependent part of epsilon and contributions C to its derivatives fac=(rrij*sigsq)**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt @@@ -737,15 -737,15 +753,15 @@@ eneps_temp(1,ij)=eneps_temp(1,ij)+e1*aux & /dabs(eps(itypi,itypj)) eneps_temp(2,ij)=eneps_temp(2,ij)+e2*aux/eps(itypi,itypj) - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij else evdw_t=evdw_t+evdwij endif if (calc_grad) then if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa write (iout,'(2(a3,i3,2x),15(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, & epsi,sigm,chi1,chi2,chip1,chip2, @@@ -792,7 -792,6 +808,7 @@@ include 'COMMON.ENEPS' include 'COMMON.IOUNITS' include 'COMMON.CALC' + include 'COMMON.SBRIDGE' logical lprn common /srutu/icall integer icant @@@ -822,6 -821,28 +838,28 @@@ C returning the ith atom to bo if (yi.lt.0) yi=yi+boxysize zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize + if ((zi.gt.bordlipbot) + &.and.(zi.lt.bordliptop)) then + C the energy transfer exist + if (zi.lt.buflipbot) then + C what fraction I am in + fracinbuf=1.0d0- + & ((zi-bordlipbot)/lipbufthick) + C lipbufthick is thickenes of lipid buffore + sslipi=sscalelip(fracinbuf) + ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick + elseif (zi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) + sslipi=sscalelip(fracinbuf) + ssgradlipi=sscagradlip(fracinbuf)/lipbufthick + else + sslipi=1.0d0 + ssgradlipi=0.0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0 + endif dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) @@@ -832,26 -853,6 +870,26 @@@ C Calculate SC interaction energy C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) + IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN + call dyn_ssbond_ene(i,j,evdwij) + evdw=evdw+evdwij +C write (iout,'(a6,2i5,0pf7.3,a3,2f10.3)') +C & 'evdw',i,j,evdwij,' ss',evdw,evdw_t +C triple bond artifac removal + do k=j+1,iend(i,iint) +C search over all next residues + if (dyn_ss_mask(k)) then +C check if they are cysteins +C write(iout,*) 'k=',k + call triple_ssbond_ene(i,j,k,evdwij) +C call the energy function that removes the artifical triple disulfide +C bond the soubroutine is located in ssMD.F + evdw=evdw+evdwij +C write (iout,'(a6,2i5,0pf7.3,a3,2f10.3)') +C & 'evdw',i,j,evdwij,'tss',evdw,evdw_t + endif!dyn_ss_mask(k) + enddo! k + ELSE ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@@ -886,6 -887,33 +924,33 @@@ C returning jth atom to bo if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize + if ((zj.gt.bordlipbot) + &.and.(zj.lt.bordliptop)) then + C the energy transfer exist + if (zj.lt.buflipbot) then + C what fraction I am in + fracinbuf=1.0d0- + & ((zj-bordlipbot)/lipbufthick) + C lipbufthick is thickenes of lipid buffore + sslipj=sscalelip(fracinbuf) + ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick + elseif (zj.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) + sslipj=sscalelip(fracinbuf) + ssgradlipj=sscagradlip(fracinbuf)/lipbufthick + else + sslipj=1.0d0 + ssgradlipj=0.0 + endif + else + sslipj=0.0d0 + ssgradlipj=0.0 + endif + aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 + & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 + bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 + C write(iout,*),aa,aa_lip(itypi,itypj),aa_aq(itypi,itypj) C checking the distance dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 xj_safe=xj @@@ -931,6 -959,7 +996,7 @@@ c write (iout,*) i,j,xj,yj,z if (sss.le.0.0) cycle C Calculate angle-dependent terms of energy and contributions to their C derivatives. + call sc_angular sigsq=1.0D0/sigsq sig=sig0ij*dsqrt(sigsq) @@@ -944,13 -973,13 +1010,13 @@@ C I hate to put IF's in the loops, but c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt evdwij=evdwij*eps2rt*eps3rt - if (bb(itypi,itypj).gt.0) then + if (bb.gt.0) then evdw=evdw+evdwij*sss else evdw_t=evdw_t+evdwij*sss @@@ -964,8 -993,8 +1030,8 @@@ c write (iout,*) "i",i," j", c & " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)), c & aux*e2/eps(itypi,itypj) c if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa #ifdef DEBUG write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, @@@ -990,8 -1019,6 +1056,8 @@@ C Calculate the radial part of the grad C Calculate angular part of the gradient. call sc_grad endif +C write(iout,*) "partial sum", evdw, evdw_t + ENDIF ! dyn_ss enddo ! j enddo ! iint enddo ! i @@@ -1097,15 -1124,15 +1163,15 @@@ C I hate to put IF's in the loops, but c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt fac_augm=rrij**expon e_augm=augm(itypi,itypj)*fac_augm evdwij=evdwij*eps2rt*eps3rt - if (bb(itypi,itypj).gt.0.0d0) then + if (bb.gt.0.0d0) then evdw=evdw+evdwij+e_augm else evdw_t=evdw_t+evdwij+e_augm @@@ -1816,6 -1843,8 +1882,8 @@@ c print *,"itilde3 i iti iti1",i do k=1,2 mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1) enddo + C write (iout,*) 'mumu',i,b1(1,iti),Ub2(1,i-2) + C Vectors and matrices dependent on a single virtual-bond dihedral. call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1)) call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2)) @@@ -2437,8 -2466,9 +2505,9 @@@ C Contribution to the local-electrostat eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) & +a33*muij(4) c write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij - c write (iout,'(a6,2i5,0pf7.3)') - c & 'eelloc',i,j,eel_loc_ij + C write (iout,'(a6,2i5,0pf7.3)') + C & 'eelloc',i,j,eel_loc_ij + C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4) c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) eel_loc=eel_loc+eel_loc_ij C Partial derivatives in virtual-bond dihedral angles gamma @@@ -2776,6 -2806,16 +2845,16 @@@ C Cartesian derivative enddo endif else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then + 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) + 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) goto 178 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Fourth-order contributions @@@ -2915,6 -2955,7 +2994,7 @@@ C Remaining derivatives of this turn co gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) enddo endif + 178 continue endif return end @@@ -3132,13 -3173,10 +3212,13 @@@ include 'COMMON.DERIV' include 'COMMON.VAR' include 'COMMON.INTERACT' + include 'COMMON.CONTROL' + include 'COMMON.IOUNITS' dimension ggg(3) ehpb=0.0D0 cd print *,'edis: nhpb=',nhpb,' fbr=',fbr cd print *,'link_start=',link_start,' link_end=',link_end +C write(iout,*) link_end, "link_end" if (link_end.eq.0) return do i=link_start,link_end C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a @@@ -3155,98 -3193,25 +3235,98 @@@ C iii and jjj point to the residues fo endif C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. - if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. +C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. +C & iabs(itype(jjj)).eq.1) then +C write(iout,*) constr_dist,"const" + if (.not.dyn_ss .and. i.le.nss) then + if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. & iabs(itype(jjj)).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij - else -C Calculate the distance between the two points and its difference from the -C target distance. - dd=dist(ii,jj) - rdis=dd-dhpb(i) + endif !ii.gt.neres + else if (ii.gt.nres .and. jj.gt.nres) then +c Restraints from contact prediction + dd=dist(ii,jj) + if (constr_dist.eq.11) then +C ehpb=ehpb+fordepth(i)**4.0d0 +C & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd +C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, +C & ehpb,fordepth(i),dd +C write(iout,*) ehpb,"atu?" +C ehpb,"tu?" +C fac=fordepth(i)**4.0d0 +C & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd + else + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "beta nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else + dd=dist(ii,jj) + rdis=dd-dhpb(i) +C Get the force constant corresponding to this distance. + waga=forcon(i) +C Calculate the contribution to energy. + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "beta reg",dd,waga*rdis*rdis +C +C Evaluate gradient. +C + fac=waga*rdis/dd + endif !end dhpb1(i).gt.0 + endif !end const_dist=11 + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo + do j=1,3 + ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) + ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) + enddo + do k=1,3 + ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) + enddo + else !ii.gt.nres +C write(iout,*) "before" + dd=dist(ii,jj) +C write(iout,*) "after",dd + if (constr_dist.eq.11) then + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd +C ehpb=ehpb+fordepth(i)**4*rlornmr1(dd,dhpb(i),dhpb1(i)) +C fac=fordepth(i)**4*rlornmr1prim(dd,dhpb(i),dhpb1(i))/dd +C print *,ehpb,"tu?" +C write(iout,*) ehpb,"btu?", +C & dd,dhpb(i),dhpb1(i),fordepth(i),forcon(i) +C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, +C & ehpb,fordepth(i),dd + else + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "alph nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else + rdis=dd-dhpb(i) C Get the force constant corresponding to this distance. - waga=forcon(i) + waga=forcon(i) C Calculate the contribution to energy. - ehpb=ehpb+waga*rdis*rdis + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "alpha reg",dd,waga*rdis*rdis C C Evaluate gradient. C - fac=waga*rdis/dd -cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, -cd & ' waga=',waga,' fac=',fac + fac=waga*rdis/dd + endif + endif + do j=1,3 ggg(j)=fac*(c(j,jj)-c(j,ii)) enddo @@@ -3266,7 -3231,7 +3346,7 @@@ C Cartesian gradient in the SC vectors enddo endif enddo - ehpb=0.5D0*ehpb + if (constr_dist.ne.11) ehpb=0.5D0*ehpb return end C-------------------------------------------------------------------------- @@@ -3455,7 -3420,7 +3535,7 @@@ c & AKSC(j,iti),abond0(j,iti), end #ifdef CRYST_THETA C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@@ -3472,14 -3437,13 +3552,14 @@@ include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.FFIELD' + include 'COMMON.TORCNSTR' common /calcthet/ term1,term2,termm,diffak,ratak, & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq, & delthe0,sig0inv,sigtc,sigsqtc,delthec,it double precision y(2),z(2) delta=0.02d0*pi - time11=dexp(-2*time) - time12=1.0d0 +c time11=dexp(-2*time) +c time12=1.0d0 etheta=0.0D0 c write (iout,*) "nres",nres c write (*,'(a,i2)') 'EBEND ICG=',icg @@@ -3512,8 -3476,8 +3592,8 @@@ C Zero the energy function and its deri if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) - icrc=0 - call proc_proc(phii,icrc) +c icrc=0 +c call proc_proc(phii,icrc) if (icrc.eq.1) phii=150.0 #else phii=phi(i) @@@ -3528,8 -3492,8 +3608,8 @@@ if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) - icrc=0 - call proc_proc(phii1,icrc) +c icrc=0 +c call proc_proc(phii1,icrc) if (icrc.eq.1) phii1=150.0 phii1=pinorm(phii1) z(1)=cos(phii1) @@@ -3599,34 -3563,7 +3679,34 @@@ c & rad2deg*phii,rad2deg*phii1,e if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett) - 1215 continue +c 1215 continue + enddo + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=1,ntheta_constr + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif +C if (energy_dec) then +C write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", +C & i,itheta,rad2deg*thetiii, +C & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), +C & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, +C & gloc(itheta+nphi-2,icg) + endif enddo C Ufff.... We've done all this!!! return @@@ -3741,7 -3678,7 +3821,7 @@@ C "Thank you" to MAPLE (probably spare end #else C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@@ -3761,7 -3698,6 +3841,7 @@@ include 'COMMON.NAMES' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.TORCNSTR' double precision coskt(mmaxtheterm),sinkt(mmaxtheterm), & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle), & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble), @@@ -3808,9 -3744,8 +3888,9 @@@ C if (itype(i-1).eq.ntyp1) cycl enddo else phii=0.0d0 - ityp1=nthetyp+1 +c ityp1=nthetyp+1 do k=1,nsingle + ityp1=ithetyp((itype(i-2))) cosph1(k)=0.0d0 sinph1(k)=0.0d0 enddo @@@ -3831,8 -3766,7 +3911,8 @@@ enddo else phii1=0.0d0 - ityp3=nthetyp+1 +c ityp3=nthetyp+1 + ityp3=ithetyp((itype(i))) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@@ -3949,36 -3883,7 +4029,36 @@@ c call flush(iout etheta=etheta+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)=wang*dethetai +c gloc(nphi+i-2,icg)=wang*dethetai + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai + enddo +C now constrains + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=1,ntheta_constr + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif +C if (energy_dec) then +C write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", +C & i,itheta,rad2deg*thetiii, +C & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), +C & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, +C & gloc(itheta+nphi-2,icg) +C endif enddo return end @@@ -4748,16 -4653,15 +4828,16 @@@ c write (iout,*) 'i=',i,' gloc=', difi=phii-phi0(i) if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 endif -! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, -! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) +C write (iout,'(a6,2i5,2f8.3,2e14.5)') "edih", +C & i,itori,rad2deg*phii, +C & rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg) enddo ! write (iout,*) 'edihcnstr',edihcnstr return @@@ -4847,24 -4751,21 +4927,24 @@@ c write (iout,*) 'i=',i,' gloc=', edihi=0.0d0 if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 - edihi=0.25d0*ftors*difi**4 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 + edihi=0.25d0*ftors(i)*difi**4 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 - edihi=0.25d0*ftors*difi**4 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 + edihi=0.25d0*ftors(i)*difi**4 else difi=0.0d0 endif + write (iout,'(a6,2i5,2f8.3,2e14.5)') "edih", + & i,itori,rad2deg*phii, + & rad2deg*difi,0.25d0*ftors(i)*difi**4 c write (iout,'(2i5,4f10.5,e15.5)') i,itori,phii,phi0(i),difi, c & drange(i),edihi ! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, -! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) +! & rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg) enddo ! write (iout,*) 'edihcnstr',edihcnstr return @@@ -5011,7 -4912,6 +5091,7 @@@ c 3 = SC...Ca...Ca...SC esccor=esccor+v1ij*cosphi+v2ij*sinphi gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) enddo +C write (iout,*)"EBACK_SC_COR",esccor,i c write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp, c & nterm_sccor(isccori,isccori1),isccori,isccori1 c gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci @@@ -7774,6 -7674,125 +7854,125 @@@ cd write (2,*) 'eel_turn6',ekont*e return end crc------------------------------------------------- + CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + subroutine Eliptransfer(eliptran) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.LOCAL' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' + include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' + include 'COMMON.SBRIDGE' + C this is done by Adasko + C print *,"wchodze" + C structure of box: + C water + C--bordliptop-- buffore starts + C--bufliptop--- here true lipid starts + C lipid + C--buflipbot--- lipid ends buffore starts + C--bordlipbot--buffore ends + eliptran=0.0 + do i=1,nres + C do i=1,1 + if (itype(i).eq.ntyp1) cycle + + positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize)) + if (positi.le.0) positi=positi+boxzsize + C print *,i + C first for peptide groups + c for each residue check if it is in lipid or lipid water border area + if ((positi.gt.bordlipbot) + &.and.(positi.lt.bordliptop)) then + C the energy transfer exist + if (positi.lt.buflipbot) then + C what fraction I am in + fracinbuf=1.0d0- + & ((positi-bordlipbot)/lipbufthick) + C lipbufthick is thickenes of lipid buffore + sslip=sscalelip(fracinbuf) + ssgradlip=-sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*pepliptran + gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0 + gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0 + C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran + elseif (positi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) + sslip=sscalelip(fracinbuf) + ssgradlip=sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*pepliptran + gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0 + gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0 + C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran + C print *, "doing sscalefor top part" + C print *,i,sslip,fracinbuf,ssgradlip + else + eliptran=eliptran+pepliptran + C print *,"I am in true lipid" + endif + C else + C eliptran=elpitran+0.0 ! I am in water + endif + enddo + C print *, "nic nie bylo w lipidzie?" + C now multiply all by the peptide group transfer factor + C eliptran=eliptran*pepliptran + C now the same for side chains + CV do i=1,1 + do i=1,nres + if (itype(i).eq.ntyp1) cycle + positi=(mod(c(3,i+nres),boxzsize)) + if (positi.le.0) positi=positi+boxzsize + C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop + c for each residue check if it is in lipid or lipid water border area + C respos=mod(c(3,i+nres),boxzsize) + C print *,positi,bordlipbot,buflipbot + if ((positi.gt.bordlipbot) + & .and.(positi.lt.bordliptop)) then + C the energy transfer exist + if (positi.lt.buflipbot) then + fracinbuf=1.0d0- + & ((positi-bordlipbot)/lipbufthick) + C lipbufthick is thickenes of lipid buffore + sslip=sscalelip(fracinbuf) + ssgradlip=-sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*liptranene(itype(i)) + gliptranx(3,i)=gliptranx(3,i) + &+ssgradlip*liptranene(itype(i)) + gliptranc(3,i-1)= gliptranc(3,i-1) + &+ssgradlip*liptranene(itype(i)) + C print *,"doing sccale for lower part" + elseif (positi.gt.bufliptop) then + fracinbuf=1.0d0- + &((bordliptop-positi)/lipbufthick) + sslip=sscalelip(fracinbuf) + ssgradlip=sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*liptranene(itype(i)) + gliptranx(3,i)=gliptranx(3,i) + &+ssgradlip*liptranene(itype(i)) + gliptranc(3,i-1)= gliptranc(3,i-1) + &+ssgradlip*liptranene(itype(i)) + C print *, "doing sscalefor top part",sslip,fracinbuf + else + eliptran=eliptran+liptranene(itype(i)) + C print *,"I am in true lipid" + endif + endif ! if in lipid or buffor + C else + C eliptran=elpitran+0.0 ! I am in water + enddo + return + end + + + CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + SUBROUTINE MATVEC2(A1,V1,V2) implicit real*8 (a-h,o-z) include 'DIMENSIONS' @@@ -7937,4 -7956,34 +8136,34 @@@ C-------------------------------------- return end C----------------------------------------------------------------------- + C----------------------------------------------------------------------- + double precision function sscalelip(r) + double precision r,gamm + include "COMMON.SPLITELE" + C if(r.lt.r_cut-rlamb) then + C sscale=1.0d0 + C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then + C gamm=(r-(r_cut-rlamb))/rlamb + sscalelip=1.0d0+r*r*(2*r-3.0d0) + C else + C sscale=0d0 + C endif + return + end + C----------------------------------------------------------------------- + double precision function sscagradlip(r) + double precision r,gamm + include "COMMON.SPLITELE" + C if(r.lt.r_cut-rlamb) then + C sscagrad=0.0d0 + C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then + C gamm=(r-(r_cut-rlamb))/rlamb + sscagradlip=r*(6*r-6.0d0) + C else + C sscagrad=0.0d0 + C endif + return + end + + C----------------------------------------------------------------------- diff --combined source/wham/src-M/initialize_p.F index c66de63,f47a951..9b03d11 --- a/source/wham/src-M/initialize_p.F +++ b/source/wham/src-M/initialize_p.F @@@ -62,6 -62,8 +62,8 @@@ ihist=30 iweight=31 izsc=32 + C Lipidic input file for parameters range 60-79 + iliptranpar=60 C C Set default weights of the energy terms. C @@@ -88,8 -90,10 +90,10 @@@ enddo do i=1,ntyp do j=1,ntyp - aa(i,j)=0.0D0 - bb(i,j)=0.0D0 + aa_lip(i,j)=0.0D0 + bb_lip(i,j)=0.0D0 + aa_aq(i,j)=0.0D0 + bb_aq(i,j)=0.0D0 augm(i,j)=0.0D0 sigma(i,j)=0.0D0 r0(i,j)=0.0D0 @@@ -184,7 -188,6 +188,7 @@@ C Initialize the bridge array do i=1,maxres ihpb(i)=0 jhpb(i)=0 + dyn_ss_mask(i)=.false. enddo C C Initialize timing. @@@ -262,19 -265,17 +266,19 @@@ c-------------------------------------- & "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ", & "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ", & "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP", - & "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T","ELT"/ + & "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T","ELIPTRAN", + & "EAFM","ETHETC"/ data wname / & "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC", & "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD", - & "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC","WLT"/ + & "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC", + & "WLIPTRAN","WAFM","WTHETC"/ data ww0 /1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0, & 1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,0.4d0,1.0d0,1.0d0, - & 0.0d0,0.0,0.0/ + & 0.0d0,0.0,0.0d0,0.0d0,0.0d0/ data nprint_ene /22/ data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19, - & 16,15,17,20,21,22/ + & 16,15,17,20,21,24,22,23/ end c--------------------------------------------------------------------------- subroutine init_int_table @@@ -331,7 -332,6 +335,7 @@@ cd write (iout,*) 'ns=',ns,' nss=',n cd & (ihpb(i),jhpb(i),i=1,nss) do i=nnt,nct-1 scheck=.false. + if (dyn_ss) goto 10 do ii=1,nss if (ihpb(ii).eq.i+nres) then scheck=.true. diff --combined source/wham/src-M/parmread.F index 0d61fee,914c090..0cf2755 --- a/source/wham/src-M/parmread.F +++ b/source/wham/src-M/parmread.F @@@ -35,8 -35,6 +35,8 @@@ character*16 key integer iparm double precision ip,mp + character*6 res1 +C write (iout,*) "KURWA" C C Body C @@@ -57,66 -55,6 +57,66 @@@ C Assign virtual-bond lengt write (iout,*) "iparm",iparm," myparm",myparm c If reading not own parameters, skip assignment + call reada(controlcard,"D0CM",d0cm,3.78d0) + call reada(controlcard,"AKCM",akcm,15.1d0) + call reada(controlcard,"AKTH",akth,11.0d0) + call reada(controlcard,"AKCT",akct,12.0d0) + call reada(controlcard,"V1SS",v1ss,-1.08d0) + call reada(controlcard,"V2SS",v2ss,7.61d0) + call reada(controlcard,"V3SS",v3ss,13.7d0) + call reada(controlcard,"EBR",ebr,-5.50D0) + call reada(controlcard,"DTRISS",dtriss,1.0D0) + call reada(controlcard,"ATRISS",atriss,0.3D0) + call reada(controlcard,"BTRISS",btriss,0.02D0) + call reada(controlcard,"CTRISS",ctriss,1.0D0) + dyn_ss=(index(controlcard,'DYN_SS').gt.0) + write(iout,*) "ATRISS",atriss + write(iout,*) "BTRISS",btriss + write(iout,*) "CTRISS",ctriss + write(iout,*) "DTRISS",dtriss + +C do i=1,maxres +C dyn_ss_mask(i)=.false. +C enddo +C ebr=-12.0D0 +c +c Old arbitrary potential - commented out. +c +c dbr= 4.20D0 +c fbr= 3.30D0 +c +c Constants of the disulfide-bond potential determined based on the RHF/6-31G** +c energy surface of diethyl disulfide. +c A. Liwo and U. Kozlowska, 11/24/03 +c + D0CM = 3.78d0 + AKCM = 15.1d0 + AKTH = 11.0d0 + AKCT = 12.0d0 + V1SS =-1.08d0 + V2SS = 7.61d0 + V3SS = 13.7d0 + + do i=1,maxres-1 + do j=i+1,maxres + dyn_ssbond_ij(i,j)=1.0d300 + enddo + enddo + call reada(controlcard,"HT",Ht,0.0D0) +C if (dyn_ss) then +C ss_depth=ebr/wsc-0.25*eps(1,1) +C write(iout,*) HT,wsc,eps(1,1),'KURWA' +C Ht=Ht/wsc-0.25*eps(1,1) + +C akcm=akcm*whpb/wsc +C akth=akth*whpb/wsc +C akct=akct*whpb/wsc +C v1ss=v1ss*whpb/wsc +C v2ss=v2ss*whpb/wsc +C v3ss=v3ss*whpb/wsc +C else +C ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb +C endif if (iparm.eq.myparm .or. .not.separate_parset) then @@@ -140,8 -78,7 +140,9 @@@ wvdwpp=ww(16) wbond=ww(18) wsccor=ww(19) + whpb=ww(15) + wstrain=ww(15) + wliptran=ww(22) endif call card_concat(controlcard,.false.) @@@ -251,6 -188,11 +252,11 @@@ enddo enddo endif + read(iliptranpar,*) pepliptran + do i=1,ntyp + read(iliptranpar,*) liptranene(i) + enddo + close(iliptranpar) #ifdef CRYST_THETA C C Read the parameters of the probability distribution/energy expression @@@ -358,7 -300,7 +364,7 @@@ C Read the parameters of Utheta determined from ab initio surfaces C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203 C -c write (iout,*) "tu dochodze" + write (iout,*) "tu dochodze" read (ithep,*) nthetyp,ntheterm,ntheterm2, & ntheterm3,nsingle,ndouble nntheterm=max0(ntheterm,ntheterm2,ntheterm3) @@@ -366,7 -308,7 +372,7 @@@ do i=-ntyp1,-1 ithetyp(i)=-ithetyp(-i) enddo -c write (iout,*) "tu dochodze" + write (iout,*) "tu dochodze" do iblock=1,2 do i=-maxthetyp,maxthetyp do j=-maxthetyp,maxthetyp @@@ -395,13 -337,11 +401,13 @@@ enddo enddo enddo +C write (iout,*) "KURWA1" do iblock=1,2 do i=0,nthetyp do j=-nthetyp,nthetyp do k=-nthetyp,nthetyp read (ithep,'(6a)') res1 + write(iout,*) res1,i,j,k read (ithep,*) aa0thet(i,j,k,iblock) read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm) read (ithep,*) @@@ -419,7 -359,6 +425,7 @@@ enddo enddo enddo +C write(iout,*) "KURWA1.1" C C For dummy ends assign glycine-type coefficients of theta-only terms; the C coefficients of theta-and-gamma-dependent terms are zero. @@@ -439,7 -378,6 +445,7 @@@ aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0 enddo enddo +C write(iout,*) "KURWA1.5" C Substitution for D aminoacids from symmetry. do iblock=1,2 do i=-nthetyp,0 @@@ -518,7 -456,7 +524,7 @@@ call flush(iout) endif #endif - +C write(iout,*) 'KURWA2' #ifdef CRYST_SC C C Read the parameters of the probability distribution/energy expression @@@ -626,7 -564,6 +632,7 @@@ enddo #endif close(irotam) +C write (iout,*) 'KURWAKURWA' #ifdef CRYST_TOR C C Read torsional parameters in old format @@@ -1101,6 -1038,13 +1107,13 @@@ C---------------------- GB or BP potent read (isidep,*)(sigii(i),i=1,ntyp) read (isidep,*)(chip(i),i=1,ntyp) read (isidep,*)(alp(i),i=1,ntyp) + do i=1,ntyp + read (isidep,*)(epslip(i,j),j=i,ntyp) + C print *,"WARNING!!" + C do j=1,ntyp + C epslip(i,j)=epslip(i,j)+0.05d0 + C enddo + enddo C For the GB potential convert sigma'**2 into chi' if (ipot.eq.4) then do i=1,ntyp @@@ -1139,6 -1083,7 +1152,7 @@@ C Calculate the "working" parameters o do i=2,ntyp do j=1,i-1 eps(i,j)=eps(j,i) + epslip(i,j)=epslip(j,i) enddo enddo do i=1,ntyp @@@ -1156,6 -1101,7 +1170,7 @@@ do i=1,ntyp do j=i,ntyp epsij=eps(i,j) + epsijlip=epslip(i,j) if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then rrij=sigma(i,j) else @@@ -1167,10 -1113,16 +1182,16 @@@ epsij=eps(i,j) sigeps=dsign(1.0D0,epsij) epsij=dabs(epsij) - aa(i,j)=epsij*rrij*rrij - bb(i,j)=-sigeps*epsij*rrij - aa(j,i)=aa(i,j) - bb(j,i)=bb(i,j) + aa_aq(i,j)=epsij*rrij*rrij + bb_aq(i,j)=-sigeps*epsij*rrij + aa_aq(j,i)=aa_aq(i,j) + bb_aq(j,i)=bb_aq(i,j) + sigeps=dsign(1.0D0,epsijlip) + epsijlip=dabs(epsijlip) + aa_lip(i,j)=epsijlip*rrij*rrij + bb_lip(i,j)=-sigeps*epsijlip*rrij + aa_lip(j,i)=aa_lip(i,j) + bb_lip(j,i)=bb_lip(i,j) if (ipot.gt.2) then sigt1sq=sigma0(i)**2 sigt2sq=sigma0(j)**2 @@@ -1259,7 -1211,7 +1280,7 @@@ C C Define the constants of the disulfide bridge C - ebr=-5.50D0 +C ebr=-12.0D0 c c Old arbitrary potential - commented out. c @@@ -1270,36 -1222,21 +1291,36 @@@ c Constants of the disulfide-bond poten c energy surface of diethyl disulfide. c A. Liwo and U. Kozlowska, 11/24/03 c - D0CM = 3.78d0 - AKCM = 15.1d0 - AKTH = 11.0d0 - AKCT = 12.0d0 - V1SS =-1.08d0 - V2SS = 7.61d0 - V3SS = 13.7d0 +C D0CM = 3.78d0 +C AKCM = 15.1d0 +C AKTH = 11.0d0 +C AKCT = 12.0d0 +C V1SS =-1.08d0 +C V2SS = 7.61d0 +C V3SS = 13.7d0 + write (iout,*) dyn_ss,'dyndyn' + if (dyn_ss) then + ss_depth=ebr/wsc-0.25*eps(1,1) +C write(iout,*) akcm,whpb,wsc,'KURWA' + Ht=Ht/wsc-0.25*eps(1,1) - if (lprint) then + akcm=akcm*whpb/wsc + akth=akth*whpb/wsc + akct=akct*whpb/wsc + v1ss=v1ss*whpb/wsc + v2ss=v2ss*whpb/wsc + v3ss=v3ss*whpb/wsc + else + ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb + endif + +C if (lprint) then write (iout,'(/a)') "Disulfide bridge parameters:" write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, & ' v3ss:',v3ss - endif +C endif return end diff --combined source/wham/src-M/readrtns.F index 557d27b,276e7a6..684f090 --- a/source/wham/src-M/readrtns.F +++ b/source/wham/src-M/readrtns.F @@@ -18,7 -18,6 +18,7 @@@ include "COMMON.CONTROL" include "COMMON.ENERGIES" include "COMMON.SPLITELE" + include "COMMON.SBRIDGE" character*800 controlcard integer i,j,k,ii,n_ene_found integer ind,itype1,itype2,itypf,itypsc,itypp @@@ -82,6 -81,23 +82,23 @@@ c Cutoff range for interactions call reada(controlcard,"R_CUT",r_cut,15.0d0) call reada(controlcard,"LAMBDA",rlamb,0.3d0) + call reada(controlcard,"LIPTHICK",lipthick,0.0d0) + call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0) + if (lipthick.gt.0.0d0) then + bordliptop=(boxzsize+lipthick)/2.0 + bordlipbot=bordliptop-lipthick + C endif + if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) + & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" + buflipbot=bordlipbot+lipbufthick + bufliptop=bordliptop-lipbufthick + if ((lipbufthick*2.0d0).gt.lipthick) + &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF" + endif + write(iout,*) "bordliptop=",bordliptop + write(iout,*) "bordlipbot=",bordlipbot + write(iout,*) "bufliptop=",bufliptop + write(iout,*) "buflipbot=",buflipbot call readi(controlcard,'SYM',symetr,1) write (iout,*) "DISTCHAINMAX",distchainmax write (iout,*) "delta",delta @@@ -98,10 -114,7 +115,10 @@@ zscfile=index(controlcard,"ZSCFILE").gt.0 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0 write (iout,*) "with_dihed_constr ",with_dihed_constr + with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0 + write (iout,*) "with_theta_constr ",with_theta_constr call readi(controlcard,'CONSTR_DIST',constr_dist,0) + dyn_ss=index(controlcard,"DYN_SS").gt.0 return end c------------------------------------------------------------------------------ @@@ -411,7 -424,7 +428,7 @@@ c-------------------------------------- external ilen,iroof double precision rmsdev,energia(0:max_ene),efree,eini,temp double precision prop(maxQ) - integer ntot_all(maxslice,0:maxprocs-1) + integer ntot_all(maxslice,0:maxprocs-1), maxslice_buff integer iparm,ib,iib,ir,nprop,nthr,npars double precision etot,time integer ixdrf,iret @@@ -542,13 -555,7 +559,13 @@@ c DA scratchfile #ifdef MPI c Check if everyone has the same number of conformations - call MPI_Allgather(stot(1),maxslice,MPI_INTEGER, + +c call MPI_ALLgather(MPI_IN_PLACE,stot(1),MPI_DATATYPE_NULL, +c & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR) + + maxslice_buff=maxslice + + call MPI_Allgather(stot(1),maxslice_buff,MPI_INTEGER, & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR) lerr=.false. do i=0,nprocs-1