source/cluster/wham/src/diff

   1 4c4
   2 <       include 'sizesclu.dat'
   3 ---
   4 >       include 'DIMENSIONS.ZSCOPT'
   5 5a6
   6 > #ifndef ISNAN
   7 6a8
   8 > #endif
   9 83a86,89
  10 > C
  11 > C 21/5/07 Calculate local sicdechain correlation energy
  12 > C
  13 >       call eback_sc_corr(esccor)
  14 99,125d104
  15 < C     call multibody(ecorr)
  16 < C
  17 < C Sum the energies
  18 < C
  19 < C scale large componenets
  20 < c#ifdef SCALE
  21 < c      ecorr5_scal=1000.0
  22 < c      eel_loc_scal=100.0
  23 < c      eello_turn3_scal=100.0
  24 < c      eello_turn4_scal=100.0
  25 < c      eturn6_scal=1000.0
  26 < c      ecorr6_scal=1000.0
  27 < c#else
  28 < c      ecorr5_scal=1.0
  29 < c      eel_loc_scal=1.0
  30 < c      eello_turn3_scal=1.0
  31 < c      eello_turn4_scal=1.0
  32 < c      eturn6_scal=1.0
  33 < c      ecorr6_scal=1.0
  34 < c#endif
  35 < c
  36 < c      ecorr5=ecorr5/ecorr5_scal
  37 < c      eel_loc=eel_loc/eel_loc_scal
  38 < c      eello_turn3=eello_turn3/eello_turn3_scal
  39 < c      eello_turn4=eello_turn4/eello_turn4_scal
  40 < c      eturn6=eturn6/eturn6_scal
  41 < c      ecorr6=ecorr6/ecorr6_scal
  42 133c112
  43 <      & +wbond*estr
  44 ---
  45 >      & +wbond*estr+wsccor*fact(1)*esccor
  46 141c120
  47 <      & +wbond*estr
  48 ---
  49 >      & +wbond*estr+wsccor*fact(1)*esccor
  50 172c151,152
  51 <       energia(19)=edihcnstr
  52 ---
  53 >       energia(19)=esccor
  54 >       energia(20)=edihcnstr
  55 173a154,160
  56 > #ifdef ISNAN
  57 > #ifdef AIX
  58 >       if (isnan(etot).ne.0) energia(0)=1.0d+99
  59 > #else
  60 >       if (isnan(etot)) energia(0)=1.0d+99
  61 > #endif
  62 > #else
  63 180a168
  64 > #endif
  65 201c189,190
  66 <      &                wturn6*fact(5)*gcorr6_turn(j,i)
  67 ---
  68 >      &                wturn6*fact(5)*gcorr6_turn(j,i)+
  69 >      &                wsccor*fact(2)*gsccorc(j,i)
  70 204c193,194
  71 <      &                  wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)
  72 ---
  73 >      &                  wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
  74 >      &                  wsccor*fact(2)*gsccorx(j,i)
  75 218c208,209
  76 <      &                wturn6*fact(5)*gcorr6_turn(j,i)
  77 ---
  78 >      &                wturn6*fact(5)*gcorr6_turn(j,i)+
  79 >      &                wsccor*fact(2)*gsccorc(j,i)
  80 221c212,213
  81 <      &                  wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)
  82 ---
  83 >      &                  wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
  84 >      &                  wsccor*fact(1)*gsccorx(j,i)
  85 224,225d215
  86 < cd      print '(i3,9(1pe12.4))',i,(gvdwc(k,i),k=1,3),(gelc(k,i),k=1,3),
  87 < cd   &        (gradc(k,i),k=1,3)
  88 230d219
  89 < cd        write (iout,*) i,g_corr5_loc(i)
  90 237a227
  91 >      &   +wsccor*fact(1)*gsccor_loc(i)
  92 240,244d229
  93 < cd    print*,evdw,wsc,evdw2,wscp,ees+evdw1,welec,ebe,wang,
  94 < cd   &  escloc,wscloc,etors,wtor,ehpb,wstrain,nss,ebr,etot
  95 < cd    call enerprint(energia(0),fact)
  96 < cd    call intout
  97 < cd    stop
  98 251c236
  99 <       include 'sizesclu.dat'
 100 ---
 101 >       include 'DIMENSIONS.ZSCOPT'
 102 279c264,265
 103 <       edihcnstr=energia(19)
 104 ---
 105 >       esccor=energia(19)
 106 >       edihcnstr=energia(20)
 107 289c275
 108 <      &  edihcnstr,ebr*nss,etot
 109 ---
 110 >      &  esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot
 111 308a295
 112 >      & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
 113 318c305,306
 114 <      &  eello_turn6,wturn6*fact(5),edihcnstr,ebr*nss,etot
 115 ---
 116 >      &  eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor,
 117 >      &  edihcnstr,ebr*nss,etot
 118 336a325
 119 >      & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
 120 351c340
 121 <       include 'sizesclu.dat'
 122 ---
 123 >       include 'DIMENSIONS.ZSCOPT'
 124 360a350
 125 >       include 'COMMON.ENEPS'
 126 368a359,363
 127 >       do i=1,210
 128 >         do j=1,2
 129 >           eneps_temp(j,i)=0.0d0
 130 >         enddo
 131 >       enddo
 132 398a394,395
 133 >             eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij)
 134 >             eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps0ij
 135 512c509
 136 <       include 'sizesclu.dat'
 137 ---
 138 >       include 'DIMENSIONS.ZSCOPT'
 139 519a517
 140 >       include 'COMMON.ENEPS'
 141 526a525,529
 142 >       do i=1,210
 143 >         do j=1,2
 144 >           eneps_temp(j,i)=0.0d0
 145 >         enddo
 146 >       enddo
 147 553a557,559
 148 >             eneps_temp(1,ij)=eneps_temp(1,ij)+(e1+a_augm)
 149 >      &        /dabs(eps(itypi,itypj))
 150 >             eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps(itypi,itypj)
 151 601c607
 152 <       include 'sizesclu.dat'
 153 ---
 154 >       include 'DIMENSIONS.ZSCOPT'
 155 609a616
 156 >       include 'COMMON.ENEPS'
 157 616a624,628
 158 >       do i=1,210
 159 >         do j=1,2
 160 >           eneps_temp(j,i)=0.0d0
 161 >         enddo
 162 >       enddo
 163 688a701,703
 164 >             eneps_temp(1,ij)=eneps_temp(1,ij)+e1*aux
 165 >      &        /dabs(eps(itypi,itypj))
 166 >             eneps_temp(2,ij)=eneps_temp(2,ij)+e2*aux/eps(itypi,itypj)
 167 728c743
 168 <       include 'sizesclu.dat'
 169 ---
 170 >       include 'DIMENSIONS.ZSCOPT'
 171 736a752
 172 >       include 'COMMON.ENEPS'
 173 742a759,763
 174 >       do i=1,210
 175 >         do j=1,2
 176 >           eneps_temp(j,i)=0.0d0
 177 >         enddo
 178 >       enddo
 179 819a841,843
 180 >             eneps_temp(1,ij)=eneps_temp(1,ij)+aux*e1
 181 >      &        /dabs(eps(itypi,itypj))
 182 >             eneps_temp(2,ij)=eneps_temp(2,ij)+aux*e2/eps(itypi,itypj)
 183 859c883
 184 <       include 'sizesclu.dat'
 185 ---
 186 >       include 'DIMENSIONS.ZSCOPT'
 187 867a892
 188 >       include 'COMMON.ENEPS'
 189 873a899,903
 190 >       do i=1,210
 191 >         do j=1,2
 192 >           eneps_temp(j,i)=0.0d0
 193 >         enddo
 194 >       enddo
 195 952a983,987
 196 >             eneps_temp(1,ij)=eneps_temp(1,ij)+aux*(e1+e_augm)
 197 >      &        /dabs(eps(itypi,itypj))
 198 >             eneps_temp(2,ij)=eneps_temp(2,ij)+aux*e2/eps(itypi,itypj)
 199 > c            eneps_temp(ij)=eneps_temp(ij)
 200 > c     &         +(evdwij+e_augm)/eps(itypi,itypj)
 201 1035c1070
 202 <       include 'sizesclu.dat'
 203 ---
 204 >       include 'DIMENSIONS.ZSCOPT'
 205 1073c1108
 206 <       include 'sizesclu.dat'
 207 ---
 208 >       include 'DIMENSIONS.ZSCOPT'
 209 1232c1267
 210 <       include 'sizesclu.dat'
 211 ---
 212 >       include 'DIMENSIONS.ZSCOPT'
 213 1415c1450
 214 <       include 'sizesclu.dat'
 215 ---
 216 >       include 'DIMENSIONS.ZSCOPT'
 217 1500c1535
 218 <       include 'sizesclu.dat'
 219 ---
 220 >       include 'DIMENSIONS.ZSCOPT'
 221 1683c1718
 222 <       include 'sizesclu.dat'
 223 ---
 224 >       include 'DIMENSIONS.ZSCOPT'
 225 2432c2467
 226 <       include 'sizesclu.dat'
 227 ---
 228 >       include 'DIMENSIONS.ZSCOPT'
 229 2699c2734
 230 <       include 'sizesclu.dat'
 231 ---
 232 >       include 'DIMENSIONS.ZSCOPT'
 233 2810c2845
 234 <       include 'sizesclu.dat'
 235 ---
 236 >       include 'DIMENSIONS.ZSCOPT'
 237 2887c2922
 238 <       include 'sizesclu.dat'
 239 ---
 240 >       include 'DIMENSIONS.ZSCOPT'
 241 2968c3003
 242 <       include 'sizesclu.dat'
 243 ---
 244 >       include 'DIMENSIONS.ZSCOPT'
 245 2978a3014
 246 >       double precision u(3),ud(3)
 247 2988,2989c3024,3027
 248 <       estr=AKP*estr
 249 < c      write (iout,*) "estr",estr
 250 ---
 251 >       estr=0.5d0*AKP*estr
 252 > c
 253 > c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
 254 > c
 255 2993,2999c3031,3070
 256 <           diff=vbld(i+nres)-vbldsc0(iti)
 257 < c          write (iout,*) i,iti,vbld(i+nres),vbldsc0(iti),diff,
 258 < c     &      AKSC(iti)*diff*diff
 259 <           estr=estr+AKSC(iti)*diff*diff
 260 <           do j=1,3
 261 <             gradbx(j,i)=AKSC(iti)*diff*dc(j,i+nres)/vbld(i+nres)
 262 <           enddo
 263 ---
 264 >           nbi=nbondterm(iti)
 265 >           if (nbi.eq.1) then
 266 >             diff=vbld(i+nres)-vbldsc0(1,iti)
 267 > c            write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
 268 > c     &      AKSC(1,iti),AKSC(1,iti)*diff*diff
 269 >             estr=estr+0.5d0*AKSC(1,iti)*diff*diff
 270 >             do j=1,3
 271 >               gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
 272 >             enddo
 273 >           else
 274 >             do j=1,nbi
 275 >               diff=vbld(i+nres)-vbldsc0(j,iti)
 276 >               ud(j)=aksc(j,iti)*diff
 277 >               u(j)=abond0(j,iti)+0.5d0*ud(j)*diff
 278 >             enddo
 279 >             uprod=u(1)
 280 >             do j=2,nbi
 281 >               uprod=uprod*u(j)
 282 >             enddo
 283 >             usum=0.0d0
 284 >             usumsqder=0.0d0
 285 >             do j=1,nbi
 286 >               uprod1=1.0d0
 287 >               uprod2=1.0d0
 288 >               do k=1,nbi
 289 >                 if (k.ne.j) then
 290 >                   uprod1=uprod1*u(k)
 291 >                   uprod2=uprod2*u(k)*u(k)
 292 >                 endif
 293 >               enddo
 294 >               usum=usum+uprod1
 295 >               usumsqder=usumsqder+ud(j)*uprod2
 296 >             enddo
 297 > c            write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti),
 298 > c     &      AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi)
 299 >             estr=estr+uprod/usum
 300 >             do j=1,3
 301 >              gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
 302 >             enddo
 303 >           endif
 304 3002,3003d3072
 305 < c      write (iout,*) "estr",estr
 306 <       estr=0.5d0*estr
 307 3005a3075
 308 > #ifdef CRYST_THETA
 309 3014c3084
 310 <       include 'sizesclu.dat'
 311 ---
 312 >       include 'DIMENSIONS.ZSCOPT'
 313 3244a3315,3509
 314 > #else
 315 > C--------------------------------------------------------------------------
 316 >       subroutine ebend(etheta)
 317 > C
 318 > C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
 319 > C angles gamma and its derivatives in consecutive thetas and gammas.
 320 > C ab initio-derived potentials from
 321 > c Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 322 > C
 323 >       implicit real*8 (a-h,o-z)
 324 >       include 'DIMENSIONS'
 325 >       include 'DIMENSIONS.ZSCOPT'
 326 >       include 'COMMON.LOCAL'
 327 >       include 'COMMON.GEO'
 328 >       include 'COMMON.INTERACT'
 329 >       include 'COMMON.DERIV'
 330 >       include 'COMMON.VAR'
 331 >       include 'COMMON.CHAIN'
 332 >       include 'COMMON.IOUNITS'
 333 >       include 'COMMON.NAMES'
 334 >       include 'COMMON.FFIELD'
 335 >       include 'COMMON.CONTROL'
 336 >       double precision coskt(mmaxtheterm),sinkt(mmaxtheterm),
 337 >      & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle),
 338 >      & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble),
 339 >      & sinph1ph2(maxdouble,maxdouble)
 340 >       logical lprn /.false./, lprn1 /.false./
 341 >       etheta=0.0D0
 342 > c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
 343 >       do i=ithet_start,ithet_end
 344 >         dethetai=0.0d0
 345 >         dephii=0.0d0
 346 >         dephii1=0.0d0
 347 >         theti2=0.5d0*theta(i)
 348 >         ityp2=ithetyp(itype(i-1))
 349 >         do k=1,nntheterm
 350 >           coskt(k)=dcos(k*theti2)
 351 >           sinkt(k)=dsin(k*theti2)
 352 >         enddo
 353 >         if (i.gt.3) then
 354 > #ifdef OSF
 355 >           phii=phi(i)
 356 >           if (phii.ne.phii) phii=150.0
 357 > #else
 358 >           phii=phi(i)
 359 > #endif
 360 >           ityp1=ithetyp(itype(i-2))
 361 >           do k=1,nsingle
 362 >             cosph1(k)=dcos(k*phii)
 363 >             sinph1(k)=dsin(k*phii)
 364 >           enddo
 365 >         else
 366 >           phii=0.0d0
 367 >           ityp1=nthetyp+1
 368 >           do k=1,nsingle
 369 >             cosph1(k)=0.0d0
 370 >             sinph1(k)=0.0d0
 371 >           enddo
 372 >         endif
 373 >         if (i.lt.nres) then
 374 > #ifdef OSF
 375 >           phii1=phi(i+1)
 376 >           if (phii1.ne.phii1) phii1=150.0
 377 >           phii1=pinorm(phii1)
 378 > #else
 379 >           phii1=phi(i+1)
 380 > #endif
 381 >           ityp3=ithetyp(itype(i))
 382 >           do k=1,nsingle
 383 >             cosph2(k)=dcos(k*phii1)
 384 >             sinph2(k)=dsin(k*phii1)
 385 >           enddo
 386 >         else
 387 >           phii1=0.0d0
 388 >           ityp3=nthetyp+1
 389 >           do k=1,nsingle
 390 >             cosph2(k)=0.0d0
 391 >             sinph2(k)=0.0d0
 392 >           enddo
 393 >         endif
 394 > c        write (iout,*) "i",i," ityp1",itype(i-2),ityp1,
 395 > c     &   " ityp2",itype(i-1),ityp2," ityp3",itype(i),ityp3
 396 > c        call flush(iout)
 397 >         ethetai=aa0thet(ityp1,ityp2,ityp3)
 398 >         do k=1,ndouble
 399 >           do l=1,k-1
 400 >             ccl=cosph1(l)*cosph2(k-l)
 401 >             ssl=sinph1(l)*sinph2(k-l)
 402 >             scl=sinph1(l)*cosph2(k-l)
 403 >             csl=cosph1(l)*sinph2(k-l)
 404 >             cosph1ph2(l,k)=ccl-ssl
 405 >             cosph1ph2(k,l)=ccl+ssl
 406 >             sinph1ph2(l,k)=scl+csl
 407 >             sinph1ph2(k,l)=scl-csl
 408 >           enddo
 409 >         enddo
 410 >         if (lprn) then
 411 >         write (iout,*) "i",i," ityp1",ityp1," ityp2",ityp2,
 412 >      &    " ityp3",ityp3," theti2",theti2," phii",phii," phii1",phii1
 413 >         write (iout,*) "coskt and sinkt"
 414 >         do k=1,nntheterm
 415 >           write (iout,*) k,coskt(k),sinkt(k)
 416 >         enddo
 417 >         endif
 418 >         do k=1,ntheterm
 419 >           ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3)*sinkt(k)
 420 >           dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3)
 421 >      &      *coskt(k)
 422 >           if (lprn)
 423 >      &    write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3),
 424 >      &     " ethetai",ethetai
 425 >         enddo
 426 >         if (lprn) then
 427 >         write (iout,*) "cosph and sinph"
 428 >         do k=1,nsingle
 429 >           write (iout,*) k,cosph1(k),sinph1(k),cosph2(k),sinph2(k)
 430 >         enddo
 431 >         write (iout,*) "cosph1ph2 and sinph2ph2"
 432 >         do k=2,ndouble
 433 >           do l=1,k-1
 434 >             write (iout,*) l,k,cosph1ph2(l,k),cosph1ph2(k,l),
 435 >      &         sinph1ph2(l,k),sinph1ph2(k,l)
 436 >           enddo
 437 >         enddo
 438 >         write(iout,*) "ethetai",ethetai
 439 >         endif
 440 >         do m=1,ntheterm2
 441 >           do k=1,nsingle
 442 >             aux=bbthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)
 443 >      &         +ccthet(k,m,ityp1,ityp2,ityp3)*sinph1(k)
 444 >      &         +ddthet(k,m,ityp1,ityp2,ityp3)*cosph2(k)
 445 >      &         +eethet(k,m,ityp1,ityp2,ityp3)*sinph2(k)
 446 >             ethetai=ethetai+sinkt(m)*aux
 447 >             dethetai=dethetai+0.5d0*m*aux*coskt(m)
 448 >             dephii=dephii+k*sinkt(m)*(
 449 >      &          ccthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)-
 450 >      &          bbthet(k,m,ityp1,ityp2,ityp3)*sinph1(k))
 451 >             dephii1=dephii1+k*sinkt(m)*(
 452 >      &          eethet(k,m,ityp1,ityp2,ityp3)*cosph2(k)-
 453 >      &          ddthet(k,m,ityp1,ityp2,ityp3)*sinph2(k))
 454 >             if (lprn)
 455 >      &      write (iout,*) "m",m," k",k," bbthet",
 456 >      &         bbthet(k,m,ityp1,ityp2,ityp3)," ccthet",
 457 >      &         ccthet(k,m,ityp1,ityp2,ityp3)," ddthet",
 458 >      &         ddthet(k,m,ityp1,ityp2,ityp3)," eethet",
 459 >      &         eethet(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai
 460 >           enddo
 461 >         enddo
 462 >         if (lprn)
 463 >      &  write(iout,*) "ethetai",ethetai
 464 >         do m=1,ntheterm3
 465 >           do k=2,ndouble
 466 >             do l=1,k-1
 467 >               aux=ffthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
 468 >      &            ffthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+
 469 >      &            ggthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
 470 >      &            ggthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)
 471 >               ethetai=ethetai+sinkt(m)*aux
 472 >               dethetai=dethetai+0.5d0*m*coskt(m)*aux
 473 >               dephii=dephii+l*sinkt(m)*(
 474 >      &           -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-
 475 >      &            ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
 476 >      &            ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
 477 >      &            ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
 478 >               dephii1=dephii1+(k-l)*sinkt(m)*(
 479 >      &           -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
 480 >      &            ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
 481 >      &            ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-
 482 >      &            ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
 483 >               if (lprn) then
 484 >               write (iout,*) "m",m," k",k," l",l," ffthet",
 485 >      &            ffthet(l,k,m,ityp1,ityp2,ityp3),
 486 >      &            ffthet(k,l,m,ityp1,ityp2,ityp3)," ggthet",
 487 >      &            ggthet(l,k,m,ityp1,ityp2,ityp3),
 488 >      &            ggthet(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai
 489 >               write (iout,*) cosph1ph2(l,k)*sinkt(m),
 490 >      &            cosph1ph2(k,l)*sinkt(m),
 491 >      &            sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
 492 >               endif
 493 >             enddo
 494 >           enddo
 495 >         enddo
 496 > 10      continue
 497 >         if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)')
 498 >      &   i,theta(i)*rad2deg,phii*rad2deg,
 499 >      &   phii1*rad2deg,ethetai
 500 >         etheta=etheta+ethetai
 501 >         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
 502 >         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
 503 >         gloc(nphi+i-2,icg)=wang*dethetai
 504 >       enddo
 505 >       return
 506 >       end
 507 > #endif
 508 > #ifdef CRYST_SC
 509 3252c3517
 510 <       include 'sizesclu.dat'
 511 ---
 512 >       include 'DIMENSIONS.ZSCOPT'
 513 3525a3791,4113
 514 > #else
 515 > c----------------------------------------------------------------------------------
 516 >       subroutine esc(escloc)
 517 > C Calculate the local energy of a side chain and its derivatives in the
 518 > C corresponding virtual-bond valence angles THETA and the spherical angles
 519 > C ALPHA and OMEGA derived from AM1 all-atom calculations.
 520 > C added by Urszula Kozlowska. 07/11/2007
 521 > C
 522 >       implicit real*8 (a-h,o-z)
 523 >       include 'DIMENSIONS'
 524 >       include 'DIMENSIONS.ZSCOPT'
 525 >       include 'COMMON.GEO'
 526 >       include 'COMMON.LOCAL'
 527 >       include 'COMMON.VAR'
 528 >       include 'COMMON.SCROT'
 529 >       include 'COMMON.INTERACT'
 530 >       include 'COMMON.DERIV'
 531 >       include 'COMMON.CHAIN'
 532 >       include 'COMMON.IOUNITS'
 533 >       include 'COMMON.NAMES'
 534 >       include 'COMMON.FFIELD'
 535 >       include 'COMMON.CONTROL'
 536 >       include 'COMMON.VECTORS'
 537 >       double precision x_prime(3),y_prime(3),z_prime(3)
 538 >      &    , sumene,dsc_i,dp2_i,x(65),
 539 >      &     xx,yy,zz,sumene1,sumene2,sumene3,sumene4,s1,s1_6,s2,s2_6,
 540 >      &    de_dxx,de_dyy,de_dzz,de_dt
 541 >       double precision s1_t,s1_6_t,s2_t,s2_6_t
 542 >       double precision
 543 >      & dXX_Ci1(3),dYY_Ci1(3),dZZ_Ci1(3),dXX_Ci(3),
 544 >      & dYY_Ci(3),dZZ_Ci(3),dXX_XYZ(3),dYY_XYZ(3),dZZ_XYZ(3),
 545 >      & dt_dCi(3),dt_dCi1(3)
 546 >       common /sccalc/ time11,time12,time112,theti,it,nlobit
 547 >       delta=0.02d0*pi
 548 >       escloc=0.0D0
 549 >       do i=loc_start,loc_end
 550 >         costtab(i+1) =dcos(theta(i+1))
 551 >         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 552 >         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 553 >         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 554 >         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 555 >         cosfac=dsqrt(cosfac2)
 556 >         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 557 >         sinfac=dsqrt(sinfac2)
 558 >         it=itype(i)
 559 >         if (it.eq.10) goto 1
 560 > c
 561 > C  Compute the axes of tghe local cartesian coordinates system; store in
 562 > c   x_prime, y_prime and z_prime
 563 > c
 564 >         do j=1,3
 565 >           x_prime(j) = 0.00
 566 >           y_prime(j) = 0.00
 567 >           z_prime(j) = 0.00
 568 >         enddo
 569 > C        write(2,*) "dc_norm", dc_norm(1,i+nres),dc_norm(2,i+nres),
 570 > C     &   dc_norm(3,i+nres)
 571 >         do j = 1,3
 572 >           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 573 >           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 574 >         enddo
 575 >         do j = 1,3
 576 >           z_prime(j) = -uz(j,i-1)
 577 >         enddo
 578 > c       write (2,*) "i",i
 579 > c       write (2,*) "x_prime",(x_prime(j),j=1,3)
 580 > c       write (2,*) "y_prime",(y_prime(j),j=1,3)
 581 > c       write (2,*) "z_prime",(z_prime(j),j=1,3)
 582 > c       write (2,*) "xx",scalar(x_prime(1),x_prime(1)),
 583 > c      & " xy",scalar(x_prime(1),y_prime(1)),
 584 > c      & " xz",scalar(x_prime(1),z_prime(1)),
 585 > c      & " yy",scalar(y_prime(1),y_prime(1)),
 586 > c      & " yz",scalar(y_prime(1),z_prime(1)),
 587 > c      & " zz",scalar(z_prime(1),z_prime(1))
 588 > c
 589 > C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 590 > C to local coordinate system. Store in xx, yy, zz.
 591 > c
 592 >         xx=0.0d0
 593 >         yy=0.0d0
 594 >         zz=0.0d0
 595 >         do j = 1,3
 596 >           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 597 >           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 598 >           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 599 >         enddo
 600 >
 601 >         xxtab(i)=xx
 602 >         yytab(i)=yy
 603 >         zztab(i)=zz
 604 > C
 605 > C Compute the energy of the ith side cbain
 606 > C
 607 > c        write (2,*) "xx",xx," yy",yy," zz",zz
 608 >         it=itype(i)
 609 >         do j = 1,65
 610 >           x(j) = sc_parmin(j,it)
 611 >         enddo
 612 > #ifdef CHECK_COORD
 613 > Cc diagnostics - remove later
 614 >         xx1 = dcos(alph(2))
 615 >         yy1 = dsin(alph(2))*dcos(omeg(2))
 616 >         zz1 = -dsin(alph(2))*dsin(omeg(2))
 617 >         write(2,'(3f8.1,3f9.3,1x,3f9.3)')
 618 >      &    alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,
 619 >      &    xx1,yy1,zz1
 620 > C,"  --- ", xx_w,yy_w,zz_w
 621 > c end diagnostics
 622 > #endif
 623 >         sumene1= x(1)+  x(2)*xx+  x(3)*yy+  x(4)*zz+  x(5)*xx**2
 624 >      &   + x(6)*yy**2+  x(7)*zz**2+  x(8)*xx*zz+  x(9)*xx*yy
 625 >      &   + x(10)*yy*zz
 626 >         sumene2=  x(11) + x(12)*xx + x(13)*yy + x(14)*zz + x(15)*xx**2
 627 >      & + x(16)*yy**2 + x(17)*zz**2 + x(18)*xx*zz + x(19)*xx*yy
 628 >      & + x(20)*yy*zz
 629 >         sumene3=  x(21) +x(22)*xx +x(23)*yy +x(24)*zz +x(25)*xx**2
 630 >      &  +x(26)*yy**2 +x(27)*zz**2 +x(28)*xx*zz +x(29)*xx*yy
 631 >      &  +x(30)*yy*zz +x(31)*xx**3 +x(32)*yy**3 +x(33)*zz**3
 632 >      &  +x(34)*(xx**2)*yy +x(35)*(xx**2)*zz +x(36)*(yy**2)*xx
 633 >      &  +x(37)*(yy**2)*zz +x(38)*(zz**2)*xx +x(39)*(zz**2)*yy
 634 >      &  +x(40)*xx*yy*zz
 635 >         sumene4= x(41) +x(42)*xx +x(43)*yy +x(44)*zz +x(45)*xx**2
 636 >      &  +x(46)*yy**2 +x(47)*zz**2 +x(48)*xx*zz +x(49)*xx*yy
 637 >      &  +x(50)*yy*zz +x(51)*xx**3 +x(52)*yy**3 +x(53)*zz**3
 638 >      &  +x(54)*(xx**2)*yy +x(55)*(xx**2)*zz +x(56)*(yy**2)*xx
 639 >      &  +x(57)*(yy**2)*zz +x(58)*(zz**2)*xx +x(59)*(zz**2)*yy
 640 >      &  +x(60)*xx*yy*zz
 641 >         dsc_i   = 0.743d0+x(61)
 642 >         dp2_i   = 1.9d0+x(62)
 643 >         dscp1=dsqrt(dsc_i**2+dp2_i**2-2*dsc_i*dp2_i
 644 >      &          *(xx*cost2tab(i+1)+yy*sint2tab(i+1)))
 645 >         dscp2=dsqrt(dsc_i**2+dp2_i**2-2*dsc_i*dp2_i
 646 >      &          *(xx*cost2tab(i+1)-yy*sint2tab(i+1)))
 647 >         s1=(1+x(63))/(0.1d0 + dscp1)
 648 >         s1_6=(1+x(64))/(0.1d0 + dscp1**6)
 649 >         s2=(1+x(65))/(0.1d0 + dscp2)
 650 >         s2_6=(1+x(65))/(0.1d0 + dscp2**6)
 651 >         sumene = ( sumene3*sint2tab(i+1) + sumene1)*(s1+s1_6)
 652 >      & + (sumene4*cost2tab(i+1) +sumene2)*(s2+s2_6)
 653 > c        write(2,'(i2," sumene",7f9.3)') i,sumene1,sumene2,sumene3,
 654 > c     &   sumene4,
 655 > c     &   dscp1,dscp2,sumene
 656 > c        sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
 657 >         escloc = escloc + sumene
 658 > c        write (2,*) "escloc",escloc
 659 >         if (.not. calc_grad) goto 1
 660 > #ifdef DEBUG
 661 > C
 662 > C This section to check the numerical derivatives of the energy of ith side
 663 > C chain in xx, yy, zz, and theta. Use the -DDEBUG compiler option or insert
 664 > C #define DEBUG in the code to turn it on.
 665 > C
 666 >         write (2,*) "sumene               =",sumene
 667 >         aincr=1.0d-7
 668 >         xxsave=xx
 669 >         xx=xx+aincr
 670 >         write (2,*) xx,yy,zz
 671 >         sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
 672 >         de_dxx_num=(sumenep-sumene)/aincr
 673 >         xx=xxsave
 674 >         write (2,*) "xx+ sumene from enesc=",sumenep
 675 >         yysave=yy
 676 >         yy=yy+aincr
 677 >         write (2,*) xx,yy,zz
 678 >         sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
 679 >         de_dyy_num=(sumenep-sumene)/aincr
 680 >         yy=yysave
 681 >         write (2,*) "yy+ sumene from enesc=",sumenep
 682 >         zzsave=zz
 683 >         zz=zz+aincr
 684 >         write (2,*) xx,yy,zz
 685 >         sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
 686 >         de_dzz_num=(sumenep-sumene)/aincr
 687 >         zz=zzsave
 688 >         write (2,*) "zz+ sumene from enesc=",sumenep
 689 >         costsave=cost2tab(i+1)
 690 >         sintsave=sint2tab(i+1)
 691 >         cost2tab(i+1)=dcos(0.5d0*(theta(i+1)+aincr))
 692 >         sint2tab(i+1)=dsin(0.5d0*(theta(i+1)+aincr))
 693 >         sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
 694 >         de_dt_num=(sumenep-sumene)/aincr
 695 >         write (2,*) " t+ sumene from enesc=",sumenep
 696 >         cost2tab(i+1)=costsave
 697 >         sint2tab(i+1)=sintsave
 698 > C End of diagnostics section.
 699 > #endif
 700 > C
 701 > C Compute the gradient of esc
 702 > C
 703 >         pom_s1=(1.0d0+x(63))/(0.1d0 + dscp1)**2
 704 >         pom_s16=6*(1.0d0+x(64))/(0.1d0 + dscp1**6)**2
 705 >         pom_s2=(1.0d0+x(65))/(0.1d0 + dscp2)**2
 706 >         pom_s26=6*(1.0d0+x(65))/(0.1d0 + dscp2**6)**2
 707 >         pom_dx=dsc_i*dp2_i*cost2tab(i+1)
 708 >         pom_dy=dsc_i*dp2_i*sint2tab(i+1)
 709 >         pom_dt1=-0.5d0*dsc_i*dp2_i*(xx*sint2tab(i+1)-yy*cost2tab(i+1))
 710 >         pom_dt2=-0.5d0*dsc_i*dp2_i*(xx*sint2tab(i+1)+yy*cost2tab(i+1))
 711 >         pom1=(sumene3*sint2tab(i+1)+sumene1)
 712 >      &     *(pom_s1/dscp1+pom_s16*dscp1**4)
 713 >         pom2=(sumene4*cost2tab(i+1)+sumene2)
 714 >      &     *(pom_s2/dscp2+pom_s26*dscp2**4)
 715 >         sumene1x=x(2)+2*x(5)*xx+x(8)*zz+ x(9)*yy
 716 >         sumene3x=x(22)+2*x(25)*xx+x(28)*zz+x(29)*yy+3*x(31)*xx**2
 717 >      &  +2*x(34)*xx*yy +2*x(35)*xx*zz +x(36)*(yy**2) +x(38)*(zz**2)
 718 >      &  +x(40)*yy*zz
 719 >         sumene2x=x(12)+2*x(15)*xx+x(18)*zz+ x(19)*yy
 720 >         sumene4x=x(42)+2*x(45)*xx +x(48)*zz +x(49)*yy +3*x(51)*xx**2
 721 >      &  +2*x(54)*xx*yy+2*x(55)*xx*zz+x(56)*(yy**2)+x(58)*(zz**2)
 722 >      &  +x(60)*yy*zz
 723 >         de_dxx =(sumene1x+sumene3x*sint2tab(i+1))*(s1+s1_6)
 724 >      &        +(sumene2x+sumene4x*cost2tab(i+1))*(s2+s2_6)
 725 >      &        +(pom1+pom2)*pom_dx
 726 > #ifdef DEBUG
 727 >         write(2,*), "de_dxx = ", de_dxx,de_dxx_num
 728 > #endif
 729 > C
 730 >         sumene1y=x(3) + 2*x(6)*yy + x(9)*xx + x(10)*zz
 731 >         sumene3y=x(23) +2*x(26)*yy +x(29)*xx +x(30)*zz +3*x(32)*yy**2
 732 >      &  +x(34)*(xx**2) +2*x(36)*yy*xx +2*x(37)*yy*zz +x(39)*(zz**2)
 733 >      &  +x(40)*xx*zz
 734 >         sumene2y=x(13) + 2*x(16)*yy + x(19)*xx + x(20)*zz
 735 >         sumene4y=x(43)+2*x(46)*yy+x(49)*xx +x(50)*zz
 736 >      &  +3*x(52)*yy**2+x(54)*xx**2+2*x(56)*yy*xx +2*x(57)*yy*zz
 737 >      &  +x(59)*zz**2 +x(60)*xx*zz
 738 >         de_dyy =(sumene1y+sumene3y*sint2tab(i+1))*(s1+s1_6)
 739 >      &        +(sumene2y+sumene4y*cost2tab(i+1))*(s2+s2_6)
 740 >      &        +(pom1-pom2)*pom_dy
 741 > #ifdef DEBUG
 742 >         write(2,*), "de_dyy = ", de_dyy,de_dyy_num
 743 > #endif
 744 > C
 745 >         de_dzz =(x(24) +2*x(27)*zz +x(28)*xx +x(30)*yy
 746 >      &  +3*x(33)*zz**2 +x(35)*xx**2 +x(37)*yy**2 +2*x(38)*zz*xx
 747 >      &  +2*x(39)*zz*yy +x(40)*xx*yy)*sint2tab(i+1)*(s1+s1_6)
 748 >      &  +(x(4) + 2*x(7)*zz+  x(8)*xx + x(10)*yy)*(s1+s1_6)
 749 >      &  +(x(44)+2*x(47)*zz +x(48)*xx   +x(50)*yy  +3*x(53)*zz**2
 750 >      &  +x(55)*xx**2 +x(57)*(yy**2)+2*x(58)*zz*xx +2*x(59)*zz*yy
 751 >      &  +x(60)*xx*yy)*cost2tab(i+1)*(s2+s2_6)
 752 >      &  + ( x(14) + 2*x(17)*zz+  x(18)*xx + x(20)*yy)*(s2+s2_6)
 753 > #ifdef DEBUG
 754 >         write(2,*), "de_dzz = ", de_dzz,de_dzz_num
 755 > #endif
 756 > C
 757 >         de_dt =  0.5d0*sumene3*cost2tab(i+1)*(s1+s1_6)
 758 >      &  -0.5d0*sumene4*sint2tab(i+1)*(s2+s2_6)
 759 >      &  +pom1*pom_dt1+pom2*pom_dt2
 760 > #ifdef DEBUG
 761 >         write(2,*), "de_dt = ", de_dt,de_dt_num
 762 > #endif
 763 > c
 764 > C
 765 >        cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres))
 766 >        cossc1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
 767 >        cosfac2xx=cosfac2*xx
 768 >        sinfac2yy=sinfac2*yy
 769 >        do k = 1,3
 770 >          dt_dCi(k) = -(dc_norm(k,i-1)+costtab(i+1)*dc_norm(k,i))*
 771 >      &      vbld_inv(i+1)
 772 >          dt_dCi1(k)= -(dc_norm(k,i)+costtab(i+1)*dc_norm(k,i-1))*
 773 >      &      vbld_inv(i)
 774 >          pom=(dC_norm(k,i+nres)-cossc*dC_norm(k,i))*vbld_inv(i+1)
 775 >          pom1=(dC_norm(k,i+nres)-cossc1*dC_norm(k,i-1))*vbld_inv(i)
 776 > c         write (iout,*) "i",i," k",k," pom",pom," pom1",pom1,
 777 > c     &    " dt_dCi",dt_dCi(k)," dt_dCi1",dt_dCi1(k)
 778 > c         write (iout,*) "dC_norm",(dC_norm(j,i),j=1,3),
 779 > c     &   (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i)
 780 >          dXX_Ci(k)=pom*cosfac-dt_dCi(k)*cosfac2xx
 781 >          dXX_Ci1(k)=-pom1*cosfac-dt_dCi1(k)*cosfac2xx
 782 >          dYY_Ci(k)=pom*sinfac+dt_dCi(k)*sinfac2yy
 783 >          dYY_Ci1(k)=pom1*sinfac+dt_dCi1(k)*sinfac2yy
 784 >          dZZ_Ci1(k)=0.0d0
 785 >          dZZ_Ci(k)=0.0d0
 786 >          do j=1,3
 787 >            dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
 788 >            dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
 789 >          enddo
 790 >
 791 >          dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
 792 >          dYY_XYZ(k)=vbld_inv(i+nres)*(y_prime(k)-yy*dC_norm(k,i+nres))
 793 >          dZZ_XYZ(k)=vbld_inv(i+nres)*(z_prime(k)-zz*dC_norm(k,i+nres))
 794 > c
 795 >          dt_dCi(k) = -dt_dCi(k)/sinttab(i+1)
 796 >          dt_dCi1(k)= -dt_dCi1(k)/sinttab(i+1)
 797 >        enddo
 798 >
 799 >        do k=1,3
 800 >          dXX_Ctab(k,i)=dXX_Ci(k)
 801 >          dXX_C1tab(k,i)=dXX_Ci1(k)
 802 >          dYY_Ctab(k,i)=dYY_Ci(k)
 803 >          dYY_C1tab(k,i)=dYY_Ci1(k)
 804 >          dZZ_Ctab(k,i)=dZZ_Ci(k)
 805 >          dZZ_C1tab(k,i)=dZZ_Ci1(k)
 806 >          dXX_XYZtab(k,i)=dXX_XYZ(k)
 807 >          dYY_XYZtab(k,i)=dYY_XYZ(k)
 808 >          dZZ_XYZtab(k,i)=dZZ_XYZ(k)
 809 >        enddo
 810 >
 811 >        do k = 1,3
 812 > c         write (iout,*) "k",k," dxx_ci1",dxx_ci1(k)," dyy_ci1",
 813 > c     &    dyy_ci1(k)," dzz_ci1",dzz_ci1(k)
 814 > c         write (iout,*) "k",k," dxx_ci",dxx_ci(k)," dyy_ci",
 815 > c     &    dyy_ci(k)," dzz_ci",dzz_ci(k)
 816 > c         write (iout,*) "k",k," dt_dci",dt_dci(k)," dt_dci",
 817 > c     &    dt_dci(k)
 818 > c         write (iout,*) "k",k," dxx_XYZ",dxx_XYZ(k)," dyy_XYZ",
 819 > c     &    dyy_XYZ(k)," dzz_XYZ",dzz_XYZ(k)
 820 >          gscloc(k,i-1)=gscloc(k,i-1)+de_dxx*dxx_ci1(k)
 821 >      &    +de_dyy*dyy_ci1(k)+de_dzz*dzz_ci1(k)+de_dt*dt_dCi1(k)
 822 >          gscloc(k,i)=gscloc(k,i)+de_dxx*dxx_Ci(k)
 823 >      &    +de_dyy*dyy_Ci(k)+de_dzz*dzz_Ci(k)+de_dt*dt_dCi(k)
 824 >          gsclocx(k,i)=                 de_dxx*dxx_XYZ(k)
 825 >      &    +de_dyy*dyy_XYZ(k)+de_dzz*dzz_XYZ(k)
 826 >        enddo
 827 > c       write(iout,*) "ENERGY GRAD = ", (gscloc(k,i-1),k=1,3),
 828 > c     &  (gscloc(k,i),k=1,3),(gsclocx(k,i),k=1,3)
 829 >
 830 > C to check gradient call subroutine check_grad
 831 >
 832 >     1 continue
 833 >       enddo
 834 >       return
 835 >       end
 836 > #endif
 837 3563c4151
 838 <       include 'sizesclu.dat'
 839 ---
 840 >       include 'DIMENSIONS.ZSCOPT'
 841 3611c4199
 842 <       include 'sizesclu.dat'
 843 ---
 844 >       include 'DIMENSIONS.ZSCOPT'
 845 3694c4282
 846 <       include 'sizesclu.dat'
 847 ---
 848 >       include 'DIMENSIONS.ZSCOPT'
 849 3780c4368
 850 <       include 'sizesclu.dat'
 851 ---
 852 >       include 'DIMENSIONS.ZSCOPT'
 853 3847a4436,4486
 854 >       subroutine eback_sc_corr(esccor)
 855 > c 7/21/2007 Correlations between the backbone-local and side-chain-local
 856 > c        conformational states; temporarily implemented as differences
 857 > c        between UNRES torsional potentials (dependent on three types of
 858 > c        residues) and the torsional potentials dependent on all 20 types
 859 > c        of residues computed from AM1 energy surfaces of terminally-blocked
 860 > c        amino-acid residues.
 861 >       implicit real*8 (a-h,o-z)
 862 >       include 'DIMENSIONS'
 863 >       include 'DIMENSIONS.ZSCOPT'
 864 >       include 'COMMON.VAR'
 865 >       include 'COMMON.GEO'
 866 >       include 'COMMON.LOCAL'
 867 >       include 'COMMON.TORSION'
 868 >       include 'COMMON.SCCOR'
 869 >       include 'COMMON.INTERACT'
 870 >       include 'COMMON.DERIV'
 871 >       include 'COMMON.CHAIN'
 872 >       include 'COMMON.NAMES'
 873 >       include 'COMMON.IOUNITS'
 874 >       include 'COMMON.FFIELD'
 875 >       include 'COMMON.CONTROL'
 876 >       logical lprn
 877 > C Set lprn=.true. for debugging
 878 >       lprn=.false.
 879 > c      lprn=.true.
 880 > c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
 881 >       esccor=0.0D0
 882 >       do i=iphi_start,iphi_end
 883 >         esccor_ii=0.0D0
 884 >         itori=itype(i-2)
 885 >         itori1=itype(i-1)
 886 >         phii=phi(i)
 887 >         gloci=0.0D0
 888 >         do j=1,nterm_sccor
 889 >           v1ij=v1sccor(j,itori,itori1)
 890 >           v2ij=v2sccor(j,itori,itori1)
 891 >           cosphi=dcos(j*phii)
 892 >           sinphi=dsin(j*phii)
 893 >           esccor=esccor+v1ij*cosphi+v2ij*sinphi
 894 >           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
 895 >         enddo
 896 >         if (lprn)
 897 >      &  write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
 898 >      &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
 899 >      &  (v1sccor(j,itori,itori1),j=1,6),(v2sccor(j,itori,itori1),j=1,6)
 900 >         gsccor_loc(i-3)=gloci
 901 >       enddo
 902 >       return
 903 >       end
 904 > c------------------------------------------------------------------------------
 905 4003c4642
 906 <       include 'sizesclu.dat'
 907 ---
 908 >       include 'DIMENSIONS.ZSCOPT'
 909 4189c4828
 910 <       include 'sizesclu.dat'
 911 ---
 912 >       include 'DIMENSIONS.ZSCOPT'
 913 4498c5137
 914 <       include 'sizesclu.dat'
 915 ---
 916 >       include 'DIMENSIONS.ZSCOPT'
 917 4565c5204
 918 <       include 'sizesclu.dat'
 919 ---
 920 >       include 'DIMENSIONS.ZSCOPT'
 921 4942c5581
 922 <       include 'sizesclu.dat'
 923 ---
 924 >       include 'DIMENSIONS.ZSCOPT'
 925 5059c5698
 926 <       include 'sizesclu.dat'
 927 ---
 928 >       include 'DIMENSIONS.ZSCOPT'
 929 5460c6099
 930 <       include 'sizesclu.dat'
 931 ---
 932 >       include 'DIMENSIONS.ZSCOPT'
 933 5597c6236
 934 <       include 'sizesclu.dat'
 935 ---
 936 >       include 'DIMENSIONS.ZSCOPT'
 937 5703c6342
 938 <       include 'sizesclu.dat'
 939 ---
 940 >       include 'DIMENSIONS.ZSCOPT'
 941 5890c6529
 942 <       include 'sizesclu.dat'
 943 ---
 944 >       include 'DIMENSIONS.ZSCOPT'
 945 6006c6645
 946 <       include 'sizesclu.dat'
 947 ---
 948 >       include 'DIMENSIONS.ZSCOPT'
 949 6252c6891
 950 <       include 'sizesclu.dat'
 951 ---
 952 >       include 'DIMENSIONS.ZSCOPT'