X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-M-SAXS-homology%2Fenergy_p_new.F;h=5360778340da83026b2b50ab92d895400debfeeb;hb=5836ecdab5a8b95f079bbf6e07374dee3fce8a26;hp=6abf7f0fed22d55dbf1326c89a6451a785582257;hpb=a05e8996b6fb955dc21aef7a85db21ee73d0a9d3;p=unres.git diff --git a/source/wham/src-M-SAXS-homology/energy_p_new.F b/source/wham/src-M-SAXS-homology/energy_p_new.F index 6abf7f0..5360778 100644 --- a/source/wham/src-M-SAXS-homology/energy_p_new.F +++ b/source/wham/src-M-SAXS-homology/energy_p_new.F @@ -124,12 +124,18 @@ C endif c print *,"Processor",myrank," computed Utord" C - call eback_sc_corr(esccor) + if (wsccor.gt.0.0d0) then + call eback_sc_corr(esccor) + else + esccor=0.0d0 + endif if (wliptran.gt.0) then call Eliptransfer(eliptran) + else + eliptran=0.0d0 endif - +#ifdef FOURBODY C C 12/1/95 Multi-body terms C @@ -151,6 +157,7 @@ c write (iout,*) ecorr,ecorr5,ecorr6,eturn6 c write (iout,*) "Calling multibody_hbond" call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1) endif +#endif c write (iout,*) "From Esaxs: Esaxs_constr",Esaxs_constr if (nsaxs.gt.0 .and. saxs_mode.eq.0) then call e_saxs(Esaxs_constr) @@ -506,10 +513,17 @@ C Bartek #ifdef SPLITELE write(iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1,wvdwpp, & estr,wbond,ebe,wang,escloc,wscloc,etors,wtor*fact(1), - & etors_d,wtor_d*fact(2),ehpb,wstrain,ecorr,wcorr*fact(3), - & ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5),eel_loc, + & etors_d,wtor_d*fact(2),ehpb,wstrain, +#ifdef FOURBODY + & ecorr,wcorr*fact(3), + & ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5), +#endif + & eel_loc, & wel_loc*fact(2),eello_turn3,wturn3*fact(2), - & eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5), + & eello_turn4,wturn4*fact(3), +#ifdef FOURBODY + & eello_turn6,wturn6*fact(5), +#endif & esccor,wsccor*fact(1),edihcnstr, & ethetacnstr,ebr*nss,Uconst,wumb,eliptran,wliptran,Eafmforc, & etube,wtube,esaxs,wsaxs,ehomology_constr, @@ -528,13 +542,17 @@ C Bartek & 'ETORSD=',1pE16.6,' WEIGHT=',1pE16.6,' (double torsional)'/ & 'EHBP= ',1pE16.6,' WEIGHT=',1pE16.6, & ' (SS bridges & dist. cnstr.)'/ +#ifdef FOURBODY & 'ECORR4=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ & 'ECORR5=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ & 'ECORR6=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ +#endif & 'EELLO= ',1pE16.6,' WEIGHT=',1pE16.6,' (electrostatic-local)'/ & 'ETURN3=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 3rd order)'/ & 'ETURN4=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 4th order)'/ +#ifdef FOURBODY & 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/ +#endif & 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (virtual-bond dihedral angle restraints)'/ & 'ETHETC=',1pE16.6,' (virtual-bond angle restraints)'/ @@ -554,10 +572,16 @@ C Bartek #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1), & estr,wbond,ebe,wang,escloc,wscloc,etors,wtor*fact(1), - & etors_d,wtor_d*fact(2),ehpb,wstrain,ecorr,wcorr*fact(3), + & etors_d,wtor_d*fact(2),ehpb,wstrain, +#ifdef FOURBODY + & ecorr,wcorr*fact(3), & ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5), +#endif & eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2), - & eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5), + & eello_turn4,wturn4*fact(3), +#ifdef FOURBODY + & eello_turn6,wturn6*fact(5), +#endif & esccor,wsccor*fact(1),edihcnstr, & ethetacnstr,ebr*nss,Uconst,wumb,eliptran,wliptran,Eafmforc, & etube,wtube,esaxs,wsaxs,ehomology_constr, @@ -575,13 +599,17 @@ C Bartek & 'ETORSD=',1pE16.6,' WEIGHT=',1pE16.6,' (double torsional)'/ & 'EHBP= ',1pE16.6,' WEIGHT=',1pE16.6, & ' (SS bridges & dist. restr.)'/ +#ifdef FOURBODY & 'ECORR4=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ & 'ECORR5=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ & 'ECORR6=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ +#endif & 'EELLO= ',1pE16.6,' WEIGHT=',1pE16.6,' (electrostatic-local)'/ & 'ETURN3=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 3rd order)'/ & 'ETURN4=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 4th order)'/ +#ifdef FOURBODY & 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/ +#endif & 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (virtual-bond dihedral angle restraints)'/ & 'ETHETC=',1pE16.6,' (virtual-bond angle restraints)'/ @@ -622,7 +650,10 @@ C include 'COMMON.SBRIDGE' include 'COMMON.NAMES' include 'COMMON.IOUNITS' +#ifdef FOURBODY include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' +#endif dimension gg(3) integer icant external icant @@ -661,6 +692,10 @@ cd & 'iend=',iend(i,iint) C Change 12/1/95 to calculate four-body interactions rij=xj*xj+yj*yj+zj*zj rrij=1.0D0/rij + sqrij=dsqrt(rij) + sss1=sscale(sqrij) + if (sss1.eq.0.0d0) cycle + sssgrad1=sscagrad(sqrij) c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj eps0ij=eps(itypi,itypj) fac=rrij**expon2 @@ -680,15 +715,16 @@ 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.gt.0.0d0) then - evdw=evdw+evdwij + evdw=evdw+sss1*evdwij else - evdw_t=evdw_t+evdwij + evdw_t=evdw_t+sss1*evdwij endif if (calc_grad) then C C Calculate the components of the gradient in DC and X C - fac=-rrij*(e1+evdwij) + fac=-rrij*(e1+evdwij)*sss1 + & +evdwij*sssgrad1/sqrij/expon gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -702,6 +738,7 @@ C enddo enddo endif +#ifdef FOURBODY C C 12/1/95, revised on 5/20/97 C @@ -758,10 +795,13 @@ cd write (iout,'(2i3,3f10.5)') cd & i,j,(gacont(kk,num_conti,i),kk=1,3) endif endif +#endif enddo ! j enddo ! iint +#ifdef FOURBODY C Change 12/1/95 num_cont(i)=num_conti +#endif enddo ! i if (calc_grad) then do i=1,nct @@ -835,6 +875,9 @@ C e_augm=augm(itypi,itypj)*fac_augm r_inv_ij=dsqrt(rrij) rij=1.0D0/r_inv_ij + sss1=sscale(rij) + if (sss1.eq.0.0d0) cycle + sssgrad1=sscagrad(rij) r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon e1=fac*fac*aa @@ -852,15 +895,16 @@ 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.gt.0.0d0) then - evdw=evdw+evdwij + evdw=evdw+evdwij*sss1 else - evdw_t=evdw_t+evdwij + evdw_t=evdw_t+evdwij*sss1 endif if (calc_grad) then C C Calculate the components of the gradient in DC and X C - fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2) + fac=(-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2))*sss1 + & +evdwij*sssgrad1*r_inv_ij/expon gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -978,6 +1022,10 @@ cd else cd rrij=rrsave(ind) cd endif rij=dsqrt(rrij) + sss1=sscale(1.0d0/rij) + if (sss1.eq.0.0d0) cycle + sssgrad1=sscagrad(1.0d0/rij) + C Calculate the angle-dependent terms of energy & contributions to derivatives. call sc_angular C Calculate whole angle-dependent part of epsilon and contributions @@ -995,9 +1043,9 @@ C to its derivatives & /dabs(eps(itypi,itypj)) eneps_temp(2,ij)=eneps_temp(2,ij)+e2*aux/eps(itypi,itypj) if (bb.gt.0.0d0) then - evdw=evdw+evdwij + evdw=evdw+sss1*evdwij else - evdw_t=evdw_t+evdwij + evdw_t=evdw_t+sss1*evdwij endif if (calc_grad) then if (lprn) then @@ -1015,6 +1063,7 @@ C Calculate gradient components. fac=-expon*(e1+evdwij) sigder=fac/sigsq fac=rrij*fac + & +evdwij*sssgrad1/sss1*rij C Calculate radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac @@ -1240,8 +1289,8 @@ C finding the closest c write (iout,*) i,j,xj,yj,zj rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) - sss=sscale((1.0d0/rij)/sigma(itypi,itypj)) - sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj)) + sss=sscale(1.0d0/rij) + sssgrad=sscagrad(1.0d0/rij) if (sss.le.0.0) cycle C Calculate angle-dependent terms of energy and contributions to their C derivatives. @@ -1401,6 +1450,9 @@ c alf12=0.0D0 dzj=dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) + sss=sscale(1.0d0/rij) + if (sss.eq.0.0d0) cycle + sssgrad=sscagrad(1.0d0/rij) C Calculate angle-dependent terms of energy and contributions to their C derivatives. call sc_angular @@ -1425,9 +1477,9 @@ c--------------------------------------------------------------- e_augm=augm(itypi,itypj)*fac_augm evdwij=evdwij*eps2rt*eps3rt if (bb.gt.0.0d0) then - evdw=evdw+evdwij+e_augm + evdw=evdw+(evdwij+e_augm)*sss else - evdw_t=evdw_t+evdwij+e_augm + evdw_t=evdw_t+(evdwij+e_augm)*sss endif ij=icant(itypi,itypj) aux=eps1*eps2rt**2*eps3rt**2 @@ -1453,6 +1505,7 @@ C Calculate gradient components. fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac-2*expon*rrij*e_augm + fac=fac+(evdwij+e_augm)*sssgrad/sss*rij C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac @@ -1730,7 +1783,7 @@ C-------------------------------------------------------------------------- include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' - include 'COMMON.CONTACTS' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' @@ -1962,6 +2015,7 @@ c & EE(1,2,iti),EE(2,2,i) c write(iout,*) "Macierz EUG", c & eug(1,1,i-2),eug(1,2,i-2),eug(2,1,i-2), c & eug(2,2,i-2) +#ifdef FOURBODY if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) & then call matmat2(CC(1,1,i-2),Ug(1,1,i-2),CUg(1,1,i-2)) @@ -1970,6 +2024,7 @@ c & eug(2,2,i-2) call matvec2(Ctilde(1,1,i-1),obrot(1,i-2),Ctobr(1,i-2)) call matvec2(Dtilde(1,1,i-2),obrot2(1,i-2),Dtobr2(1,i-2)) endif +#endif else do k=1,2 Ub2(k,i-2)=0.0d0 @@ -2011,6 +2066,7 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then #endif cd write (iout,*) 'mu1',mu1(:,i-2) cd write (iout,*) 'mu2',mu2(:,i-2) +#ifdef FOURBODY if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) & then if (calc_grad) then @@ -2033,7 +2089,9 @@ C Vectors and matrices dependent on a single virtual-bond dihedral. call matmat2(EUgder(1,1,i-2),DD(1,1,i-1),EUgDder(1,1,i-2)) endif endif +#endif enddo +#ifdef FOURBODY C Matrices dependent on two consecutive virtual-bond dihedrals. C The order of matrices is from left to right. if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) @@ -2053,6 +2111,7 @@ C The order of matrices is from left to right. endif enddo endif +#endif return end C-------------------------------------------------------------------------- @@ -2078,7 +2137,11 @@ C include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' +#ifdef FOURBODY include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' +#endif + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' @@ -2151,9 +2214,11 @@ cd enddo eello_turn3=0.0d0 eello_turn4=0.0d0 ind=0 +#ifdef FOURBODY do i=1,nres num_cont_hb(i)=0 enddo +#endif cd print '(a)','Enter EELEC' cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e do i=1,nres @@ -2204,7 +2269,9 @@ c end if num_conti=0 call eelecij(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) +#ifdef FOURBODY num_cont_hb(i)=num_conti +#endif enddo do i=iturn4_start,iturn4_end if (i.lt.1) cycle @@ -2259,13 +2326,16 @@ c endif if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize - +#ifdef FOURBODY num_conti=num_cont_hb(i) +#endif c write(iout,*) "JESTEM W PETLI" call eelecij(i,i+3,ees,evdw1,eel_loc) if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) & call eturn4(i,eello_turn4) +#ifdef FOURBODY num_cont_hb(i)=num_conti +#endif enddo ! i C Loop over all neighbouring boxes C do xshift=-1,1 @@ -2332,7 +2402,9 @@ c go to 166 c endif c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) +#ifdef FOURBODY num_conti=num_cont_hb(i) +#endif C I TU KURWA do j=ielstart(i),ielend(i) C do j=16,17 @@ -2348,7 +2420,9 @@ c & .or.itype(j-1).eq.ntyp1 &) cycle call eelecij(i,j,ees,evdw1,eel_loc) enddo ! j +#ifdef FOURBODY num_cont_hb(i)=num_conti +#endif enddo ! i C enddo ! zshift C enddo ! yshift @@ -2380,7 +2454,11 @@ C------------------------------------------------------------------------------- include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.INTERACT' +#ifdef FOURBODY include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' +#endif + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VECTORS' include 'COMMON.FFIELD' @@ -2498,8 +2576,9 @@ C yj=yj-ymedi C zj=zj-zmedi rij=xj*xj+yj*yj+zj*zj - sss=sscale(sqrt(rij)) - sssgrad=sscagrad(sqrt(rij)) + sss=sscale(sqrt(rij)) + if (sss.eq.0.0d0) return + sssgrad=sscagrad(sqrt(rij)) c write (iout,*) "ij",i,j," rij",sqrt(rij)," r_cut",r_cut, c & " rlamb",rlamb," sss",sss c if (sss.gt.0.0d0) then @@ -2667,9 +2746,10 @@ cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo if (sss.gt.0.0) then - ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj - ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj - ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj + facvdw=facvdw+sssgrad*rmij*evdwij + ggg(1)=facvdw*xj + ggg(2)=facvdw*yj + ggg(3)=facvdw*zj else ggg(1)=0.0 ggg(2)=0.0 @@ -2696,10 +2776,11 @@ cgrad enddo endif ! calc_grad #else C MARYSIA - facvdw=(ev1+evdwij)*sss + facvdw=(ev1+evdwij) facel=(el1+eesij) fac1=fac - fac=-3*rrmij*(facvdw+facvdw+facel) + fac=-3*rrmij*(facvdw+facvdw+facel)*sss + & +(evdwij+eesij)*sssgrad*rrmij erij(1)=xj*rmij erij(2)=yj*rmij erij(3)=zj*rmij @@ -3015,7 +3096,7 @@ C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eel_loc_ij=eel_loc_ij - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij c if (eel_loc_ij.ne.0) @@ -3079,7 +3160,7 @@ C Calculate patrial derivative for theta angle & +a23*gmuij1(2) & +a32*gmuij1(3) & +a33*gmuij1(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss c write(iout,*) "derivative over thatai" c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), c & a33*gmuij1(4) @@ -3095,7 +3176,7 @@ c & a33*gmuij2(4) & +a33*gmuij2(4) gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss c Derivative over j residue geel_loc_ji=a22*gmuji1(1) @@ -3120,7 +3201,7 @@ c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), c & a33*gmuji2(4) gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss #endif cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij @@ -3136,10 +3217,14 @@ C Partial derivatives in virtual-bond dihedral angles gamma & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)) & *fac_shield(i)*fac_shield(j) C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) + aux=eel_loc_ij/sss*sssgrad*rmij + ggg(1)=aux*xj + ggg(2)=aux*yj + ggg(3)=aux*zj do l=1,3 - ggg(l)=(agg(l,1)*muij(1)+ + ggg(l)=ggg(l)+(agg(l,1)*muij(1)+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l) gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l) cgrad ghalf=0.5d0*ggg(l) @@ -3176,6 +3261,7 @@ C Remaining derivatives of eello C Change 12/26/95 to calculate four-body contributions to H-bonding energy c if (j.gt.i+1 .and. num_conti.le.maxconts) then +#ifdef FOURBODY if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0 & .and. num_conti.le.maxconts) then c write (iout,*) i,j," entered corr" @@ -3315,11 +3401,17 @@ cd fprimcont=0.0D0 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k) enddo gggp(1)=gggp(1)+ees0pijp*xj + & +ees0p(num_conti,i)/sss*rmij*xj*sssgrad gggp(2)=gggp(2)+ees0pijp*yj + & +ees0p(num_conti,i)/sss*rmij*yj*sssgrad gggp(3)=gggp(3)+ees0pijp*zj + & +ees0p(num_conti,i)/sss*rmij*zj*sssgrad gggm(1)=gggm(1)+ees0mijp*xj + & +ees0m(num_conti,i)/sss*rmij*xj*sssgrad gggm(2)=gggm(2)+ees0mijp*yj + & +ees0m(num_conti,i)/sss*rmij*yj*sssgrad gggm(3)=gggm(3)+ees0mijp*zj + & +ees0m(num_conti,i)/sss*rmij*zj*sssgrad C Derivatives due to the contact function gacont_hbr(1,num_conti,i)=fprimcont*xj gacont_hbr(2,num_conti,i)=fprimcont*yj @@ -3334,28 +3426,28 @@ cgrad ghalfm=0.5D0*gggm(k) gacontp_hb1(k,num_conti,i)=!ghalfp & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gacontp_hb2(k,num_conti,i)=!ghalfp & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gacontp_hb3(k,num_conti,i)=gggp(k) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gacontm_hb1(k,num_conti,i)=!ghalfm & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gacontm_hb2(k,num_conti,i)=!ghalfm & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gacontm_hb3(k,num_conti,i)=gggm(k) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss enddo C Diagnostics. Comment out or remove after debugging! @@ -3374,6 +3466,7 @@ cdiag enddo endif ! num_conti.le.maxconts endif ! fcont.gt.0 endif ! j.gt.i+1 +#endif if (calc_grad) then if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then do k=1,4 @@ -3415,6 +3508,7 @@ C Third- and fourth-order contributions from turns include 'COMMON.FFIELD' include 'COMMON.CONTROL' include 'COMMON.SHIELD' + include 'COMMON.CORRMAT' dimension ggg(3) double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2), & e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2), @@ -3603,6 +3697,7 @@ C Third- and fourth-order contributions from turns include 'COMMON.FFIELD' include 'COMMON.CONTROL' include 'COMMON.SHIELD' + include 'COMMON.CORRMAT' dimension ggg(3) double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2), & e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2), @@ -6934,6 +7029,7 @@ c gsccor_loc(i-3)=gloci enddo return end +#ifdef FOURBODY c------------------------------------------------------------------------------ subroutine multibody(ecorr) C This subroutine calculates multi-body contributions to energy following @@ -6946,6 +7042,8 @@ C contribution equal to sqrt(eps(i,j)*eps(i+1,j+1)) is added. include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' double precision gx(3),gx1(3) logical lprn @@ -7000,6 +7098,8 @@ c------------------------------------------------------------------------------ include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' double precision gx(3),gx1(3) logical lprn lprn=.false. @@ -7042,6 +7142,8 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' double precision gx(3),gx1(3) logical lprn,ldone @@ -7115,6 +7217,8 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.CHAIN' include 'COMMON.CONTROL' include 'COMMON.SHIELD' @@ -7272,6 +7376,8 @@ c------------------------------------------------------------------------------ include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.SHIELD' include 'COMMON.CONTROL' double precision gx(3),gx1(3) @@ -7448,6 +7554,8 @@ C--------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -7514,6 +7622,8 @@ C include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -7893,6 +8003,8 @@ C--------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -8008,6 +8120,8 @@ C--------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -8425,6 +8539,8 @@ c-------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -8568,6 +8684,8 @@ c-------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -8675,6 +8793,8 @@ c---------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -8863,6 +8983,8 @@ c---------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -8981,6 +9103,8 @@ c---------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -9228,6 +9352,8 @@ c---------------------------------------------------------------------------- include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.CONTMAT' + include 'COMMON.CORRMAT' include 'COMMON.TORSION' include 'COMMON.VAR' include 'COMMON.GEO' @@ -9548,7 +9674,7 @@ cd write (2,*) 'ekont',ekont cd write (2,*) 'eel_turn6',ekont*eel_turn6 return end - +#endif crc------------------------------------------------- CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine Eliptransfer(eliptran)