+++ /dev/null
-4c4
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-5a6
-> #ifndef ISNAN
-6a8
-> #endif
-83a86,89
-> C
-> C 21/5/07 Calculate local sicdechain correlation energy
-> C
-> call eback_sc_corr(esccor)
-99,125d104
-< C call multibody(ecorr)
-< C
-< C Sum the energies
-< C
-< C scale large componenets
-< c#ifdef SCALE
-< c ecorr5_scal=1000.0
-< c eel_loc_scal=100.0
-< c eello_turn3_scal=100.0
-< c eello_turn4_scal=100.0
-< c eturn6_scal=1000.0
-< c ecorr6_scal=1000.0
-< c#else
-< c ecorr5_scal=1.0
-< c eel_loc_scal=1.0
-< c eello_turn3_scal=1.0
-< c eello_turn4_scal=1.0
-< c eturn6_scal=1.0
-< c ecorr6_scal=1.0
-< c#endif
-< c
-< c ecorr5=ecorr5/ecorr5_scal
-< c eel_loc=eel_loc/eel_loc_scal
-< c eello_turn3=eello_turn3/eello_turn3_scal
-< c eello_turn4=eello_turn4/eello_turn4_scal
-< c eturn6=eturn6/eturn6_scal
-< c ecorr6=ecorr6/ecorr6_scal
-133c112
-< & +wbond*estr
----
-> & +wbond*estr+wsccor*fact(1)*esccor
-141c120
-< & +wbond*estr
----
-> & +wbond*estr+wsccor*fact(1)*esccor
-172c151,152
-< energia(19)=edihcnstr
----
-> energia(19)=esccor
-> energia(20)=edihcnstr
-173a154,160
-> #ifdef ISNAN
-> #ifdef AIX
-> if (isnan(etot).ne.0) energia(0)=1.0d+99
-> #else
-> if (isnan(etot)) energia(0)=1.0d+99
-> #endif
-> #else
-180a168
-> #endif
-201c189,190
-< & wturn6*fact(5)*gcorr6_turn(j,i)
----
-> & wturn6*fact(5)*gcorr6_turn(j,i)+
-> & wsccor*fact(2)*gsccorc(j,i)
-204c193,194
-< & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)
----
-> & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
-> & wsccor*fact(2)*gsccorx(j,i)
-218c208,209
-< & wturn6*fact(5)*gcorr6_turn(j,i)
----
-> & wturn6*fact(5)*gcorr6_turn(j,i)+
-> & wsccor*fact(2)*gsccorc(j,i)
-221c212,213
-< & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)
----
-> & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
-> & wsccor*fact(1)*gsccorx(j,i)
-224,225d215
-< cd print '(i3,9(1pe12.4))',i,(gvdwc(k,i),k=1,3),(gelc(k,i),k=1,3),
-< cd & (gradc(k,i),k=1,3)
-230d219
-< cd write (iout,*) i,g_corr5_loc(i)
-237a227
-> & +wsccor*fact(1)*gsccor_loc(i)
-240,244d229
-< cd print*,evdw,wsc,evdw2,wscp,ees+evdw1,welec,ebe,wang,
-< cd & escloc,wscloc,etors,wtor,ehpb,wstrain,nss,ebr,etot
-< cd call enerprint(energia(0),fact)
-< cd call intout
-< cd stop
-251c236
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-279c264,265
-< edihcnstr=energia(19)
----
-> esccor=energia(19)
-> edihcnstr=energia(20)
-289c275
-< & edihcnstr,ebr*nss,etot
----
-> & esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot
-308a295
-> & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
-318c305,306
-< & eello_turn6,wturn6*fact(5),edihcnstr,ebr*nss,etot
----
-> & eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor,
-> & edihcnstr,ebr*nss,etot
-336a325
-> & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
-351c340
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-360a350
-> include 'COMMON.ENEPS'
-368a359,363
-> do i=1,210
-> do j=1,2
-> eneps_temp(j,i)=0.0d0
-> enddo
-> enddo
-398a394,395
-> eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij)
-> eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps0ij
-512c509
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-519a517
-> include 'COMMON.ENEPS'
-526a525,529
-> do i=1,210
-> do j=1,2
-> eneps_temp(j,i)=0.0d0
-> enddo
-> enddo
-553a557,559
-> eneps_temp(1,ij)=eneps_temp(1,ij)+(e1+a_augm)
-> & /dabs(eps(itypi,itypj))
-> eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps(itypi,itypj)
-601c607
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-609a616
-> include 'COMMON.ENEPS'
-616a624,628
-> do i=1,210
-> do j=1,2
-> eneps_temp(j,i)=0.0d0
-> enddo
-> enddo
-688a701,703
-> 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)
-728c743
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-736a752
-> include 'COMMON.ENEPS'
-742a759,763
-> do i=1,210
-> do j=1,2
-> eneps_temp(j,i)=0.0d0
-> enddo
-> enddo
-819a841,843
-> eneps_temp(1,ij)=eneps_temp(1,ij)+aux*e1
-> & /dabs(eps(itypi,itypj))
-> eneps_temp(2,ij)=eneps_temp(2,ij)+aux*e2/eps(itypi,itypj)
-859c883
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-867a892
-> include 'COMMON.ENEPS'
-873a899,903
-> do i=1,210
-> do j=1,2
-> eneps_temp(j,i)=0.0d0
-> enddo
-> enddo
-952a983,987
-> eneps_temp(1,ij)=eneps_temp(1,ij)+aux*(e1+e_augm)
-> & /dabs(eps(itypi,itypj))
-> eneps_temp(2,ij)=eneps_temp(2,ij)+aux*e2/eps(itypi,itypj)
-> c eneps_temp(ij)=eneps_temp(ij)
-> c & +(evdwij+e_augm)/eps(itypi,itypj)
-1035c1070
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-1073c1108
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-1232c1267
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-1415c1450
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-1500c1535
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-1683c1718
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-2432c2467
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-2699c2734
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-2810c2845
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-2887c2922
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-2968c3003
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-2978a3014
-> double precision u(3),ud(3)
-2988,2989c3024,3027
-< estr=AKP*estr
-< c write (iout,*) "estr",estr
----
-> estr=0.5d0*AKP*estr
-> c
-> c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
-> c
-2993,2999c3031,3070
-< diff=vbld(i+nres)-vbldsc0(iti)
-< c write (iout,*) i,iti,vbld(i+nres),vbldsc0(iti),diff,
-< c & AKSC(iti)*diff*diff
-< estr=estr+AKSC(iti)*diff*diff
-< do j=1,3
-< gradbx(j,i)=AKSC(iti)*diff*dc(j,i+nres)/vbld(i+nres)
-< enddo
----
-> nbi=nbondterm(iti)
-> if (nbi.eq.1) then
-> diff=vbld(i+nres)-vbldsc0(1,iti)
-> c write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
-> c & AKSC(1,iti),AKSC(1,iti)*diff*diff
-> estr=estr+0.5d0*AKSC(1,iti)*diff*diff
-> do j=1,3
-> gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
-> enddo
-> else
-> do j=1,nbi
-> diff=vbld(i+nres)-vbldsc0(j,iti)
-> ud(j)=aksc(j,iti)*diff
-> u(j)=abond0(j,iti)+0.5d0*ud(j)*diff
-> enddo
-> uprod=u(1)
-> do j=2,nbi
-> uprod=uprod*u(j)
-> enddo
-> usum=0.0d0
-> usumsqder=0.0d0
-> do j=1,nbi
-> uprod1=1.0d0
-> uprod2=1.0d0
-> do k=1,nbi
-> if (k.ne.j) then
-> uprod1=uprod1*u(k)
-> uprod2=uprod2*u(k)*u(k)
-> endif
-> enddo
-> usum=usum+uprod1
-> usumsqder=usumsqder+ud(j)*uprod2
-> enddo
-> c write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti),
-> c & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi)
-> estr=estr+uprod/usum
-> do j=1,3
-> gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
-> enddo
-> endif
-3002,3003d3072
-< c write (iout,*) "estr",estr
-< estr=0.5d0*estr
-3005a3075
-> #ifdef CRYST_THETA
-3014c3084
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-3244a3315,3509
-> #else
-> C--------------------------------------------------------------------------
-> subroutine ebend(etheta)
-> C
-> C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
-> C angles gamma and its derivatives in consecutive thetas and gammas.
-> C ab initio-derived potentials from
-> c Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
-> C
-> implicit real*8 (a-h,o-z)
-> include 'DIMENSIONS'
-> include 'DIMENSIONS.ZSCOPT'
-> include 'COMMON.LOCAL'
-> include 'COMMON.GEO'
-> include 'COMMON.INTERACT'
-> include 'COMMON.DERIV'
-> include 'COMMON.VAR'
-> include 'COMMON.CHAIN'
-> include 'COMMON.IOUNITS'
-> include 'COMMON.NAMES'
-> include 'COMMON.FFIELD'
-> include 'COMMON.CONTROL'
-> double precision coskt(mmaxtheterm),sinkt(mmaxtheterm),
-> & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle),
-> & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble),
-> & sinph1ph2(maxdouble,maxdouble)
-> logical lprn /.false./, lprn1 /.false./
-> etheta=0.0D0
-> c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
-> do i=ithet_start,ithet_end
-> dethetai=0.0d0
-> dephii=0.0d0
-> dephii1=0.0d0
-> theti2=0.5d0*theta(i)
-> ityp2=ithetyp(itype(i-1))
-> do k=1,nntheterm
-> coskt(k)=dcos(k*theti2)
-> sinkt(k)=dsin(k*theti2)
-> enddo
-> if (i.gt.3) then
-> #ifdef OSF
-> phii=phi(i)
-> if (phii.ne.phii) phii=150.0
-> #else
-> phii=phi(i)
-> #endif
-> ityp1=ithetyp(itype(i-2))
-> do k=1,nsingle
-> cosph1(k)=dcos(k*phii)
-> sinph1(k)=dsin(k*phii)
-> enddo
-> else
-> phii=0.0d0
-> ityp1=nthetyp+1
-> do k=1,nsingle
-> cosph1(k)=0.0d0
-> sinph1(k)=0.0d0
-> enddo
-> endif
-> if (i.lt.nres) then
-> #ifdef OSF
-> phii1=phi(i+1)
-> if (phii1.ne.phii1) phii1=150.0
-> phii1=pinorm(phii1)
-> #else
-> phii1=phi(i+1)
-> #endif
-> ityp3=ithetyp(itype(i))
-> do k=1,nsingle
-> cosph2(k)=dcos(k*phii1)
-> sinph2(k)=dsin(k*phii1)
-> enddo
-> else
-> phii1=0.0d0
-> ityp3=nthetyp+1
-> do k=1,nsingle
-> cosph2(k)=0.0d0
-> sinph2(k)=0.0d0
-> enddo
-> endif
-> c write (iout,*) "i",i," ityp1",itype(i-2),ityp1,
-> c & " ityp2",itype(i-1),ityp2," ityp3",itype(i),ityp3
-> c call flush(iout)
-> ethetai=aa0thet(ityp1,ityp2,ityp3)
-> do k=1,ndouble
-> do l=1,k-1
-> ccl=cosph1(l)*cosph2(k-l)
-> ssl=sinph1(l)*sinph2(k-l)
-> scl=sinph1(l)*cosph2(k-l)
-> csl=cosph1(l)*sinph2(k-l)
-> cosph1ph2(l,k)=ccl-ssl
-> cosph1ph2(k,l)=ccl+ssl
-> sinph1ph2(l,k)=scl+csl
-> sinph1ph2(k,l)=scl-csl
-> enddo
-> enddo
-> if (lprn) then
-> write (iout,*) "i",i," ityp1",ityp1," ityp2",ityp2,
-> & " ityp3",ityp3," theti2",theti2," phii",phii," phii1",phii1
-> write (iout,*) "coskt and sinkt"
-> do k=1,nntheterm
-> write (iout,*) k,coskt(k),sinkt(k)
-> enddo
-> endif
-> do k=1,ntheterm
-> ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3)*sinkt(k)
-> dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3)
-> & *coskt(k)
-> if (lprn)
-> & write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3),
-> & " ethetai",ethetai
-> enddo
-> if (lprn) then
-> write (iout,*) "cosph and sinph"
-> do k=1,nsingle
-> write (iout,*) k,cosph1(k),sinph1(k),cosph2(k),sinph2(k)
-> enddo
-> write (iout,*) "cosph1ph2 and sinph2ph2"
-> do k=2,ndouble
-> do l=1,k-1
-> write (iout,*) l,k,cosph1ph2(l,k),cosph1ph2(k,l),
-> & sinph1ph2(l,k),sinph1ph2(k,l)
-> enddo
-> enddo
-> write(iout,*) "ethetai",ethetai
-> endif
-> do m=1,ntheterm2
-> do k=1,nsingle
-> aux=bbthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)
-> & +ccthet(k,m,ityp1,ityp2,ityp3)*sinph1(k)
-> & +ddthet(k,m,ityp1,ityp2,ityp3)*cosph2(k)
-> & +eethet(k,m,ityp1,ityp2,ityp3)*sinph2(k)
-> ethetai=ethetai+sinkt(m)*aux
-> dethetai=dethetai+0.5d0*m*aux*coskt(m)
-> dephii=dephii+k*sinkt(m)*(
-> & ccthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)-
-> & bbthet(k,m,ityp1,ityp2,ityp3)*sinph1(k))
-> dephii1=dephii1+k*sinkt(m)*(
-> & eethet(k,m,ityp1,ityp2,ityp3)*cosph2(k)-
-> & ddthet(k,m,ityp1,ityp2,ityp3)*sinph2(k))
-> if (lprn)
-> & write (iout,*) "m",m," k",k," bbthet",
-> & bbthet(k,m,ityp1,ityp2,ityp3)," ccthet",
-> & ccthet(k,m,ityp1,ityp2,ityp3)," ddthet",
-> & ddthet(k,m,ityp1,ityp2,ityp3)," eethet",
-> & eethet(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai
-> enddo
-> enddo
-> if (lprn)
-> & write(iout,*) "ethetai",ethetai
-> do m=1,ntheterm3
-> do k=2,ndouble
-> do l=1,k-1
-> aux=ffthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
-> & ffthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+
-> & ggthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
-> & ggthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)
-> ethetai=ethetai+sinkt(m)*aux
-> dethetai=dethetai+0.5d0*m*coskt(m)*aux
-> dephii=dephii+l*sinkt(m)*(
-> & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-
-> & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
-> & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
-> & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
-> dephii1=dephii1+(k-l)*sinkt(m)*(
-> & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
-> & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
-> & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-
-> & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
-> if (lprn) then
-> write (iout,*) "m",m," k",k," l",l," ffthet",
-> & ffthet(l,k,m,ityp1,ityp2,ityp3),
-> & ffthet(k,l,m,ityp1,ityp2,ityp3)," ggthet",
-> & ggthet(l,k,m,ityp1,ityp2,ityp3),
-> & ggthet(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai
-> write (iout,*) cosph1ph2(l,k)*sinkt(m),
-> & cosph1ph2(k,l)*sinkt(m),
-> & sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
-> endif
-> enddo
-> enddo
-> enddo
-> 10 continue
-> if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)')
-> & i,theta(i)*rad2deg,phii*rad2deg,
-> & phii1*rad2deg,ethetai
-> 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
-> enddo
-> return
-> end
-> #endif
-> #ifdef CRYST_SC
-3252c3517
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-3525a3791,4113
-> #else
-> c----------------------------------------------------------------------------------
-> subroutine esc(escloc)
-> C Calculate the local energy of a side chain and its derivatives in the
-> C corresponding virtual-bond valence angles THETA and the spherical angles
-> C ALPHA and OMEGA derived from AM1 all-atom calculations.
-> C added by Urszula Kozlowska. 07/11/2007
-> C
-> implicit real*8 (a-h,o-z)
-> include 'DIMENSIONS'
-> include 'DIMENSIONS.ZSCOPT'
-> include 'COMMON.GEO'
-> include 'COMMON.LOCAL'
-> include 'COMMON.VAR'
-> include 'COMMON.SCROT'
-> include 'COMMON.INTERACT'
-> include 'COMMON.DERIV'
-> include 'COMMON.CHAIN'
-> include 'COMMON.IOUNITS'
-> include 'COMMON.NAMES'
-> include 'COMMON.FFIELD'
-> include 'COMMON.CONTROL'
-> include 'COMMON.VECTORS'
-> double precision x_prime(3),y_prime(3),z_prime(3)
-> & , sumene,dsc_i,dp2_i,x(65),
-> & xx,yy,zz,sumene1,sumene2,sumene3,sumene4,s1,s1_6,s2,s2_6,
-> & de_dxx,de_dyy,de_dzz,de_dt
-> double precision s1_t,s1_6_t,s2_t,s2_6_t
-> double precision
-> & dXX_Ci1(3),dYY_Ci1(3),dZZ_Ci1(3),dXX_Ci(3),
-> & dYY_Ci(3),dZZ_Ci(3),dXX_XYZ(3),dYY_XYZ(3),dZZ_XYZ(3),
-> & dt_dCi(3),dt_dCi1(3)
-> common /sccalc/ time11,time12,time112,theti,it,nlobit
-> delta=0.02d0*pi
-> escloc=0.0D0
-> do i=loc_start,loc_end
-> costtab(i+1) =dcos(theta(i+1))
-> sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
-> cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
-> sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
-> cosfac2=0.5d0/(1.0d0+costtab(i+1))
-> cosfac=dsqrt(cosfac2)
-> sinfac2=0.5d0/(1.0d0-costtab(i+1))
-> sinfac=dsqrt(sinfac2)
-> it=itype(i)
-> if (it.eq.10) goto 1
-> c
-> C Compute the axes of tghe local cartesian coordinates system; store in
-> c x_prime, y_prime and z_prime
-> c
-> do j=1,3
-> x_prime(j) = 0.00
-> y_prime(j) = 0.00
-> z_prime(j) = 0.00
-> enddo
-> C write(2,*) "dc_norm", dc_norm(1,i+nres),dc_norm(2,i+nres),
-> C & dc_norm(3,i+nres)
-> do j = 1,3
-> x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
-> y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
-> enddo
-> do j = 1,3
-> z_prime(j) = -uz(j,i-1)
-> enddo
-> c write (2,*) "i",i
-> c write (2,*) "x_prime",(x_prime(j),j=1,3)
-> c write (2,*) "y_prime",(y_prime(j),j=1,3)
-> c write (2,*) "z_prime",(z_prime(j),j=1,3)
-> c write (2,*) "xx",scalar(x_prime(1),x_prime(1)),
-> c & " xy",scalar(x_prime(1),y_prime(1)),
-> c & " xz",scalar(x_prime(1),z_prime(1)),
-> c & " yy",scalar(y_prime(1),y_prime(1)),
-> c & " yz",scalar(y_prime(1),z_prime(1)),
-> c & " zz",scalar(z_prime(1),z_prime(1))
-> c
-> C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
-> C to local coordinate system. Store in xx, yy, zz.
-> c
-> xx=0.0d0
-> yy=0.0d0
-> zz=0.0d0
-> do j = 1,3
-> xx = xx + x_prime(j)*dc_norm(j,i+nres)
-> yy = yy + y_prime(j)*dc_norm(j,i+nres)
-> zz = zz + z_prime(j)*dc_norm(j,i+nres)
-> enddo
->
-> xxtab(i)=xx
-> yytab(i)=yy
-> zztab(i)=zz
-> C
-> C Compute the energy of the ith side cbain
-> C
-> c write (2,*) "xx",xx," yy",yy," zz",zz
-> it=itype(i)
-> do j = 1,65
-> x(j) = sc_parmin(j,it)
-> enddo
-> #ifdef CHECK_COORD
-> Cc diagnostics - remove later
-> xx1 = dcos(alph(2))
-> yy1 = dsin(alph(2))*dcos(omeg(2))
-> zz1 = -dsin(alph(2))*dsin(omeg(2))
-> write(2,'(3f8.1,3f9.3,1x,3f9.3)')
-> & alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,
-> & xx1,yy1,zz1
-> C," --- ", xx_w,yy_w,zz_w
-> c end diagnostics
-> #endif
-> sumene1= x(1)+ x(2)*xx+ x(3)*yy+ x(4)*zz+ x(5)*xx**2
-> & + x(6)*yy**2+ x(7)*zz**2+ x(8)*xx*zz+ x(9)*xx*yy
-> & + x(10)*yy*zz
-> sumene2= x(11) + x(12)*xx + x(13)*yy + x(14)*zz + x(15)*xx**2
-> & + x(16)*yy**2 + x(17)*zz**2 + x(18)*xx*zz + x(19)*xx*yy
-> & + x(20)*yy*zz
-> sumene3= x(21) +x(22)*xx +x(23)*yy +x(24)*zz +x(25)*xx**2
-> & +x(26)*yy**2 +x(27)*zz**2 +x(28)*xx*zz +x(29)*xx*yy
-> & +x(30)*yy*zz +x(31)*xx**3 +x(32)*yy**3 +x(33)*zz**3
-> & +x(34)*(xx**2)*yy +x(35)*(xx**2)*zz +x(36)*(yy**2)*xx
-> & +x(37)*(yy**2)*zz +x(38)*(zz**2)*xx +x(39)*(zz**2)*yy
-> & +x(40)*xx*yy*zz
-> sumene4= x(41) +x(42)*xx +x(43)*yy +x(44)*zz +x(45)*xx**2
-> & +x(46)*yy**2 +x(47)*zz**2 +x(48)*xx*zz +x(49)*xx*yy
-> & +x(50)*yy*zz +x(51)*xx**3 +x(52)*yy**3 +x(53)*zz**3
-> & +x(54)*(xx**2)*yy +x(55)*(xx**2)*zz +x(56)*(yy**2)*xx
-> & +x(57)*(yy**2)*zz +x(58)*(zz**2)*xx +x(59)*(zz**2)*yy
-> & +x(60)*xx*yy*zz
-> dsc_i = 0.743d0+x(61)
-> dp2_i = 1.9d0+x(62)
-> dscp1=dsqrt(dsc_i**2+dp2_i**2-2*dsc_i*dp2_i
-> & *(xx*cost2tab(i+1)+yy*sint2tab(i+1)))
-> dscp2=dsqrt(dsc_i**2+dp2_i**2-2*dsc_i*dp2_i
-> & *(xx*cost2tab(i+1)-yy*sint2tab(i+1)))
-> s1=(1+x(63))/(0.1d0 + dscp1)
-> s1_6=(1+x(64))/(0.1d0 + dscp1**6)
-> s2=(1+x(65))/(0.1d0 + dscp2)
-> s2_6=(1+x(65))/(0.1d0 + dscp2**6)
-> sumene = ( sumene3*sint2tab(i+1) + sumene1)*(s1+s1_6)
-> & + (sumene4*cost2tab(i+1) +sumene2)*(s2+s2_6)
-> c write(2,'(i2," sumene",7f9.3)') i,sumene1,sumene2,sumene3,
-> c & sumene4,
-> c & dscp1,dscp2,sumene
-> c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
-> escloc = escloc + sumene
-> c write (2,*) "escloc",escloc
-> if (.not. calc_grad) goto 1
-> #ifdef DEBUG
-> C
-> C This section to check the numerical derivatives of the energy of ith side
-> C chain in xx, yy, zz, and theta. Use the -DDEBUG compiler option or insert
-> C #define DEBUG in the code to turn it on.
-> C
-> write (2,*) "sumene =",sumene
-> aincr=1.0d-7
-> xxsave=xx
-> xx=xx+aincr
-> write (2,*) xx,yy,zz
-> sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
-> de_dxx_num=(sumenep-sumene)/aincr
-> xx=xxsave
-> write (2,*) "xx+ sumene from enesc=",sumenep
-> yysave=yy
-> yy=yy+aincr
-> write (2,*) xx,yy,zz
-> sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
-> de_dyy_num=(sumenep-sumene)/aincr
-> yy=yysave
-> write (2,*) "yy+ sumene from enesc=",sumenep
-> zzsave=zz
-> zz=zz+aincr
-> write (2,*) xx,yy,zz
-> sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
-> de_dzz_num=(sumenep-sumene)/aincr
-> zz=zzsave
-> write (2,*) "zz+ sumene from enesc=",sumenep
-> costsave=cost2tab(i+1)
-> sintsave=sint2tab(i+1)
-> cost2tab(i+1)=dcos(0.5d0*(theta(i+1)+aincr))
-> sint2tab(i+1)=dsin(0.5d0*(theta(i+1)+aincr))
-> sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
-> de_dt_num=(sumenep-sumene)/aincr
-> write (2,*) " t+ sumene from enesc=",sumenep
-> cost2tab(i+1)=costsave
-> sint2tab(i+1)=sintsave
-> C End of diagnostics section.
-> #endif
-> C
-> C Compute the gradient of esc
-> C
-> pom_s1=(1.0d0+x(63))/(0.1d0 + dscp1)**2
-> pom_s16=6*(1.0d0+x(64))/(0.1d0 + dscp1**6)**2
-> pom_s2=(1.0d0+x(65))/(0.1d0 + dscp2)**2
-> pom_s26=6*(1.0d0+x(65))/(0.1d0 + dscp2**6)**2
-> pom_dx=dsc_i*dp2_i*cost2tab(i+1)
-> pom_dy=dsc_i*dp2_i*sint2tab(i+1)
-> pom_dt1=-0.5d0*dsc_i*dp2_i*(xx*sint2tab(i+1)-yy*cost2tab(i+1))
-> pom_dt2=-0.5d0*dsc_i*dp2_i*(xx*sint2tab(i+1)+yy*cost2tab(i+1))
-> pom1=(sumene3*sint2tab(i+1)+sumene1)
-> & *(pom_s1/dscp1+pom_s16*dscp1**4)
-> pom2=(sumene4*cost2tab(i+1)+sumene2)
-> & *(pom_s2/dscp2+pom_s26*dscp2**4)
-> sumene1x=x(2)+2*x(5)*xx+x(8)*zz+ x(9)*yy
-> sumene3x=x(22)+2*x(25)*xx+x(28)*zz+x(29)*yy+3*x(31)*xx**2
-> & +2*x(34)*xx*yy +2*x(35)*xx*zz +x(36)*(yy**2) +x(38)*(zz**2)
-> & +x(40)*yy*zz
-> sumene2x=x(12)+2*x(15)*xx+x(18)*zz+ x(19)*yy
-> sumene4x=x(42)+2*x(45)*xx +x(48)*zz +x(49)*yy +3*x(51)*xx**2
-> & +2*x(54)*xx*yy+2*x(55)*xx*zz+x(56)*(yy**2)+x(58)*(zz**2)
-> & +x(60)*yy*zz
-> de_dxx =(sumene1x+sumene3x*sint2tab(i+1))*(s1+s1_6)
-> & +(sumene2x+sumene4x*cost2tab(i+1))*(s2+s2_6)
-> & +(pom1+pom2)*pom_dx
-> #ifdef DEBUG
-> write(2,*), "de_dxx = ", de_dxx,de_dxx_num
-> #endif
-> C
-> sumene1y=x(3) + 2*x(6)*yy + x(9)*xx + x(10)*zz
-> sumene3y=x(23) +2*x(26)*yy +x(29)*xx +x(30)*zz +3*x(32)*yy**2
-> & +x(34)*(xx**2) +2*x(36)*yy*xx +2*x(37)*yy*zz +x(39)*(zz**2)
-> & +x(40)*xx*zz
-> sumene2y=x(13) + 2*x(16)*yy + x(19)*xx + x(20)*zz
-> sumene4y=x(43)+2*x(46)*yy+x(49)*xx +x(50)*zz
-> & +3*x(52)*yy**2+x(54)*xx**2+2*x(56)*yy*xx +2*x(57)*yy*zz
-> & +x(59)*zz**2 +x(60)*xx*zz
-> de_dyy =(sumene1y+sumene3y*sint2tab(i+1))*(s1+s1_6)
-> & +(sumene2y+sumene4y*cost2tab(i+1))*(s2+s2_6)
-> & +(pom1-pom2)*pom_dy
-> #ifdef DEBUG
-> write(2,*), "de_dyy = ", de_dyy,de_dyy_num
-> #endif
-> C
-> de_dzz =(x(24) +2*x(27)*zz +x(28)*xx +x(30)*yy
-> & +3*x(33)*zz**2 +x(35)*xx**2 +x(37)*yy**2 +2*x(38)*zz*xx
-> & +2*x(39)*zz*yy +x(40)*xx*yy)*sint2tab(i+1)*(s1+s1_6)
-> & +(x(4) + 2*x(7)*zz+ x(8)*xx + x(10)*yy)*(s1+s1_6)
-> & +(x(44)+2*x(47)*zz +x(48)*xx +x(50)*yy +3*x(53)*zz**2
-> & +x(55)*xx**2 +x(57)*(yy**2)+2*x(58)*zz*xx +2*x(59)*zz*yy
-> & +x(60)*xx*yy)*cost2tab(i+1)*(s2+s2_6)
-> & + ( x(14) + 2*x(17)*zz+ x(18)*xx + x(20)*yy)*(s2+s2_6)
-> #ifdef DEBUG
-> write(2,*), "de_dzz = ", de_dzz,de_dzz_num
-> #endif
-> C
-> de_dt = 0.5d0*sumene3*cost2tab(i+1)*(s1+s1_6)
-> & -0.5d0*sumene4*sint2tab(i+1)*(s2+s2_6)
-> & +pom1*pom_dt1+pom2*pom_dt2
-> #ifdef DEBUG
-> write(2,*), "de_dt = ", de_dt,de_dt_num
-> #endif
-> c
-> C
-> cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres))
-> cossc1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
-> cosfac2xx=cosfac2*xx
-> sinfac2yy=sinfac2*yy
-> do k = 1,3
-> dt_dCi(k) = -(dc_norm(k,i-1)+costtab(i+1)*dc_norm(k,i))*
-> & vbld_inv(i+1)
-> dt_dCi1(k)= -(dc_norm(k,i)+costtab(i+1)*dc_norm(k,i-1))*
-> & vbld_inv(i)
-> pom=(dC_norm(k,i+nres)-cossc*dC_norm(k,i))*vbld_inv(i+1)
-> pom1=(dC_norm(k,i+nres)-cossc1*dC_norm(k,i-1))*vbld_inv(i)
-> c write (iout,*) "i",i," k",k," pom",pom," pom1",pom1,
-> c & " dt_dCi",dt_dCi(k)," dt_dCi1",dt_dCi1(k)
-> c write (iout,*) "dC_norm",(dC_norm(j,i),j=1,3),
-> c & (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i)
-> dXX_Ci(k)=pom*cosfac-dt_dCi(k)*cosfac2xx
-> dXX_Ci1(k)=-pom1*cosfac-dt_dCi1(k)*cosfac2xx
-> dYY_Ci(k)=pom*sinfac+dt_dCi(k)*sinfac2yy
-> dYY_Ci1(k)=pom1*sinfac+dt_dCi1(k)*sinfac2yy
-> dZZ_Ci1(k)=0.0d0
-> dZZ_Ci(k)=0.0d0
-> do j=1,3
-> dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
-> dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
-> enddo
->
-> dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
-> dYY_XYZ(k)=vbld_inv(i+nres)*(y_prime(k)-yy*dC_norm(k,i+nres))
-> dZZ_XYZ(k)=vbld_inv(i+nres)*(z_prime(k)-zz*dC_norm(k,i+nres))
-> c
-> dt_dCi(k) = -dt_dCi(k)/sinttab(i+1)
-> dt_dCi1(k)= -dt_dCi1(k)/sinttab(i+1)
-> enddo
->
-> do k=1,3
-> dXX_Ctab(k,i)=dXX_Ci(k)
-> dXX_C1tab(k,i)=dXX_Ci1(k)
-> dYY_Ctab(k,i)=dYY_Ci(k)
-> dYY_C1tab(k,i)=dYY_Ci1(k)
-> dZZ_Ctab(k,i)=dZZ_Ci(k)
-> dZZ_C1tab(k,i)=dZZ_Ci1(k)
-> dXX_XYZtab(k,i)=dXX_XYZ(k)
-> dYY_XYZtab(k,i)=dYY_XYZ(k)
-> dZZ_XYZtab(k,i)=dZZ_XYZ(k)
-> enddo
->
-> do k = 1,3
-> c write (iout,*) "k",k," dxx_ci1",dxx_ci1(k)," dyy_ci1",
-> c & dyy_ci1(k)," dzz_ci1",dzz_ci1(k)
-> c write (iout,*) "k",k," dxx_ci",dxx_ci(k)," dyy_ci",
-> c & dyy_ci(k)," dzz_ci",dzz_ci(k)
-> c write (iout,*) "k",k," dt_dci",dt_dci(k)," dt_dci",
-> c & dt_dci(k)
-> c write (iout,*) "k",k," dxx_XYZ",dxx_XYZ(k)," dyy_XYZ",
-> c & dyy_XYZ(k)," dzz_XYZ",dzz_XYZ(k)
-> gscloc(k,i-1)=gscloc(k,i-1)+de_dxx*dxx_ci1(k)
-> & +de_dyy*dyy_ci1(k)+de_dzz*dzz_ci1(k)+de_dt*dt_dCi1(k)
-> gscloc(k,i)=gscloc(k,i)+de_dxx*dxx_Ci(k)
-> & +de_dyy*dyy_Ci(k)+de_dzz*dzz_Ci(k)+de_dt*dt_dCi(k)
-> gsclocx(k,i)= de_dxx*dxx_XYZ(k)
-> & +de_dyy*dyy_XYZ(k)+de_dzz*dzz_XYZ(k)
-> enddo
-> c write(iout,*) "ENERGY GRAD = ", (gscloc(k,i-1),k=1,3),
-> c & (gscloc(k,i),k=1,3),(gsclocx(k,i),k=1,3)
->
-> C to check gradient call subroutine check_grad
->
-> 1 continue
-> enddo
-> return
-> end
-> #endif
-3563c4151
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-3611c4199
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-3694c4282
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-3780c4368
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-3847a4436,4486
-> subroutine eback_sc_corr(esccor)
-> c 7/21/2007 Correlations between the backbone-local and side-chain-local
-> c conformational states; temporarily implemented as differences
-> c between UNRES torsional potentials (dependent on three types of
-> c residues) and the torsional potentials dependent on all 20 types
-> c of residues computed from AM1 energy surfaces of terminally-blocked
-> c amino-acid residues.
-> implicit real*8 (a-h,o-z)
-> include 'DIMENSIONS'
-> include 'DIMENSIONS.ZSCOPT'
-> include 'COMMON.VAR'
-> include 'COMMON.GEO'
-> include 'COMMON.LOCAL'
-> include 'COMMON.TORSION'
-> include 'COMMON.SCCOR'
-> include 'COMMON.INTERACT'
-> include 'COMMON.DERIV'
-> include 'COMMON.CHAIN'
-> include 'COMMON.NAMES'
-> include 'COMMON.IOUNITS'
-> include 'COMMON.FFIELD'
-> include 'COMMON.CONTROL'
-> logical lprn
-> C Set lprn=.true. for debugging
-> lprn=.false.
-> c lprn=.true.
-> c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
-> esccor=0.0D0
-> do i=iphi_start,iphi_end
-> esccor_ii=0.0D0
-> itori=itype(i-2)
-> itori1=itype(i-1)
-> phii=phi(i)
-> gloci=0.0D0
-> do j=1,nterm_sccor
-> v1ij=v1sccor(j,itori,itori1)
-> v2ij=v2sccor(j,itori,itori1)
-> cosphi=dcos(j*phii)
-> sinphi=dsin(j*phii)
-> esccor=esccor+v1ij*cosphi+v2ij*sinphi
-> gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
-> enddo
-> if (lprn)
-> & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
-> & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
-> & (v1sccor(j,itori,itori1),j=1,6),(v2sccor(j,itori,itori1),j=1,6)
-> gsccor_loc(i-3)=gloci
-> enddo
-> return
-> end
-> c------------------------------------------------------------------------------
-4003c4642
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-4189c4828
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-4498c5137
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-4565c5204
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-4942c5581
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-5059c5698
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-5460c6099
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-5597c6236
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-5703c6342
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-5890c6529
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-6006c6645
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
-6252c6891
-< include 'sizesclu.dat'
----
-> include 'DIMENSIONS.ZSCOPT'
projects / unres.git / blobdiff