source/maxlik/src_MD_T_maxlik-NEWCORR.safe/w2x_eps.F

   1       subroutine w2x(nvarr,x,*)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5       include "COMMON.CONTROL"
   6       include "COMMON.INTERACT"
   7       include "COMMON.TORSION"
   8       include "COMMON.WEIGHTS"
   9       include "COMMON.ENERGIES"
  10       include "COMMON.IOUNITS"
  11       include "COMMON.NAMES"
  12       include "COMMON.VAR"
  13       include "COMMON.VMCPAR"
  14       include "COMMON.CLASSES"
  15       include "COMMON.WEIGHTDER"
  16       include "COMMON.FFIELD"
  17       include "COMMON.SCCOR"
  18       integer itor2typ(-2:2) /-20,-20,10,9,20/
  19
  20       double precision x(max_paropt),xchuj,epsi
  21       integer nvarr,i,ii,j,k,l,ll,iprot,iblock
  22       integer iti,itj
  23
  24       ii=0
  25       do i=1,n_ene
  26         if (imask(i).eq.1) then
  27           ii=ii+1
  28           x(ii)=ww(i)
  29           x_low(ii)=ww_low(i)
  30           x_up(ii)=ww_up(i)
  31         endif
  32       enddo
  33 c
  34 c SC parameters
  35 c Change 9/21/07 AL Variables are epsilons and not the increments on epsilons
  36 c to enable control on their deviations from original values.
  37 c Single-body epsilons are now defined as epsilons of pairs of same type side
  38 c chains.
  39 c
  40       write (iout,*) "W2X: mod_side",mod_side," nsingle_sc",nsingle_sc,
  41      &  " npair_sc",npair_sc
  42       call flush(iout)
  43       if (mod_side) then
  44         do i=1,ntyp
  45           do j=1,ntyp
  46             eps_orig(i,j)=eps(i,j)
  47           enddo
  48         enddo
  49         do i=1,nsingle_sc
  50           iti=ityp_ssc(i)
  51           ii=ii+1
  52           epsi=eps(iti,iti)
  53           x(ii)=epsi
  54 c          write (iout,*) "ii",ii," x",x(ii)
  55           x_low(ii)=epss_low(i)
  56           x_up(i)=epss_up(i)
  57 c          if (epss_up(i).eq.0.0d0) then
  58 c            x_up(ii)=dabs(epsi)*epss_low(i)
  59 c            x_low(ii)=epsi-x_up(ii)
  60 c            x_up(ii)=epsi+x_up(ii)
  61 c          else
  62 c            x_low(ii)=epss_low(i)
  63 c            x_up(ii)=epss_up(i)
  64 c          endif
  65         enddo
  66         do i=1,npair_sc
  67           ii=ii+1
  68           iti=ityp_psc(1,i)
  69           itj=ityp_psc(2,i)
  70           if (iti.eq.itj) then
  71             do j=1,nsingle_sc
  72               if (ityp_ssc(j).eq.iti) then
  73                 write (iout,*) "Error - pair ",restyp(iti),"-",
  74      &           restyp(iti),
  75      &           " specified in variable epsilon list but type defined",
  76      &           " in single-type list."
  77                 return1
  78               endif
  79             enddo
  80           endif
  81           epsi=eps(iti,itj)
  82           x(ii)=epsi
  83           x_low(ii)=epsp_low(i)
  84           x_up(ii)=epsp_up(i)
  85 c          if (epsp_up(i).eq.0.0d0) then
  86 c            x_up(ii)=dabs(epsi)*epsp_low(i)
  87 c            x_low(ii)=epsi-x_up(ii)
  88 c            x_up(ii)=epsi+x_up(ii)
  89 c          else
  90 c            x_low(ii)=epsp_low(i)
  91 c            x_up(ii)=epsp_up(i)
  92 c          endif
  93 c          write (iout,*) "ii",ii," x",x(ii)
  94         enddo
  95       endif
  96
  97
  98       ntor_var=0
  99       nbacktor_var=0
 100
 101       if (mod_tor) then
 102
 103       write (iout,*) "Indices of optimizable torsionals"
 104       do i=-ntortyp+1,ntortyp-1
 105         do j=-ntortyp+1,ntortyp-1
 106 c          write (iout,*) "mask",i,j,mask_tor(j,i)
 107
 108           if (tor_mode.eq.0) then
 109
 110             do iblock=1,2
 111             if (mask_tor(0,j,i,iblock).eq.1) then
 112               ii=ii+1
 113               ntor_var=ntor_var+1
 114               x(ii)=weitor(0,j,i,iblock)
 115               x_low(ii)=weitor_low(0,j,i,iblock)
 116               x_up(ii)=weitor_up(0,j,i,iblock)
 117               write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
 118      &        restyp(itor2typ(i)),iblock,ii
 119             else if (mask_tor(0,j,i,iblock).eq.2) then
 120               do k=1,nterm(j,i,iblock)
 121                 ii=ii+1
 122                 ntor_var=ntor_var+1
 123                 x(ii)=v1(k,j,i,iblock)
 124                 x_low(ii)=weitor_low(0,j,i,iblock)
 125                 x_up(ii)=weitor_up(0,j,i,iblock)
 126                 write (iout,'(2a4,3i4)') restyp(itor2typ(j)),
 127      &          restyp(itor2typ(i)),iblock,k,ii
 128               enddo
 129             else if (mask_tor(0,j,i,iblock).gt.2) then
 130               do k=1,nterm(j,i,iblock)
 131                 ii=ii+1
 132                 ntor_var=ntor_var+1
 133                 x(ii)=v1(k,j,i,iblock)
 134                 x_low(ii)=weitor_low(0,j,i,iblock)
 135                 x_up(ii)=weitor_up(0,j,i,iblock)
 136                 write (iout,'(2a4,3i4)') restyp(itor2typ(j)),
 137      &          restyp(itor2typ(i)),iblock,k,ii
 138               enddo
 139               do k=1,nterm(j,i,iblock)
 140                 ii=ii+1
 141                 ntor_var=ntor_var+1
 142                 x(ii)=v2(k,j,i,iblock)
 143                 x_low(ii)=weitor_low(0,j,i,iblock)
 144                 x_up(ii)=weitor_up(0,j,i,iblock)
 145                 write (iout,'(2a4,3i4)') restyp(itor2typ(j)),
 146      &          restyp(itor2typ(i)),ii,k,ii
 147               enddo
 148             endif
 149             enddo ! iblock
 150
 151           else if (tor_mode.eq.1) then
 152
 153             if (mask_tor(0,j,i,1).eq.1) then
 154               ii=ii+1
 155               ntor_var=ntor_var+1
 156               x(ii)=weitor(0,j,i,1)
 157               x_low(ii)=weitor_low(0,j,i,1)
 158               x_up(ii)=weitor_up(0,j,i,1)
 159               write (iout,'(2a4,i4)') restyp(itor2typ(j)),
 160      &          restyp(itor2typ(i)),ii
 161             else if (mask_tor(0,j,i,1).eq.2) then
 162               do k=1,nterm_kcc(j,i)
 163                 do l=1,nterm_kcc_Tb(j,i)
 164                   do ll=1,nterm_kcc_Tb(j,i)
 165                     ii=ii+1
 166                     ntor_var=ntor_var+1
 167                     x(ii)=v1_kcc(ll,l,k,j,i)
 168                     x_low(ii)=weitor_low(0,j,i,1)
 169                     x_up(ii)=weitor_up(0,j,i,1)
 170                     write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
 171      &                restyp(itor2typ(i)),k,ii
 172                   enddo ! ll
 173                 enddo ! l
 174               enddo
 175             else if (mask_tor(0,j,i,1).gt.2) then
 176               do k=1,nterm_kcc(j,i)
 177                 do l=1,nterm_kcc_Tb(j,i)
 178                   do ll=1,nterm_kcc_Tb(j,i)
 179                     ii=ii+1
 180                     ntor_var=ntor_var+1
 181                     x(ii)=v1_kcc(ll,l,k,j,i)
 182                     x_low(ii)=weitor_low(0,j,i,1)
 183                     x_up(ii)=weitor_up(0,j,i,1)
 184                     write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
 185      &              restyp(itor2typ(i)),k,ii
 186                   enddo ! ll
 187                 enddo ! l
 188               enddo
 189               do k=1,nterm_kcc(j,i)
 190                 do l=1,nterm_kcc_Tb(j,i)
 191                   do ll=1,nterm_kcc_Tb(j,i)
 192                     ii=ii+1
 193                     ntor_var=ntor_var+1
 194                     x(ii)=v2_kcc(ll,l,k,j,i)
 195                     x_low(ii)=weitor_low(0,j,i,1)
 196                     x_up(ii)=weitor_up(0,j,i,1)
 197                     write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
 198      &              restyp(itor2typ(i)),k,ii
 199                   enddo! ll
 200                 enddo ! l
 201               enddo
 202             endif
 203           else
 204 #ifdef NEWCORR
 205 c Compute torsional coefficients from local energy surface expansion coeffs
 206           if (mask_tor(0,j,i,1).eq.1) then
 207               ii=ii+1
 208               write (iout,*) "i",i," j",j," weitor",weitor(0,j,i,1)
 209               x(ii)=weitor(0,j,i,1)
 210               x_low(ii)=weitor_low(0,j,i,1)
 211               x_up(ii)=weitor_up(0,j,i,1)
 212           endif
 213 #else
 214             write (iout,*) "ERROR: TOR_MODE>1 ALLOWED ONLY WITH NEWCORR"
 215             stop
 216 #endif
 217           endif ! tor_mode
 218         enddo
 219       enddo
 220
 221       endif
 222
 223       nbacktor_var=ntor_var
 224       write (iout,*) "nbacktor_var",nbacktor_var
 225
 226       if (mod_sccor) then
 227
 228       write (iout,*) "Indices of optimizable SCCOR torsionals"
 229       do i=-nsccortyp,nsccortyp
 230         do j=-nsccortyp,nsccortyp
 231           do k=1,3
 232 c          write (iout,*) "mask",i,j,mask_tor(j,i)
 233           if (mask_tor(k,j,i,1).eq.1) then
 234             ii=ii+1
 235             ntor_var=ntor_var+1
 236             x(ii)=weitor(k,j,i,1)
 237             x_low(ii)=weitor_low(k,j,i,1)
 238             x_up(ii)=weitor_up(k,j,i,1)
 239             write (iout,'(i4,2a4,i4)') k,restyp(j),restyp(i),ii
 240           else if (mask_tor(k,j,i,1).eq.2) then
 241             do l=1,nterm_sccor(j,i)
 242               ii=ii+1
 243               ntor_var=ntor_var+1
 244               x(ii)=v1sccor(l,k,j,i)
 245               x_low(ii)=weitor_low(k,j,i,1)
 246               x_up(ii)=weitor_up(k,j,i,1)
 247               write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii
 248             enddo
 249           else if (mask_tor(k,j,i,1).gt.2) then
 250             do l=1,nterm_sccor(j,i)
 251               ii=ii+1
 252               ntor_var=ntor_var+1
 253               x(ii)=v1sccor(l,k,j,i)
 254               x_low(ii)=weitor_low(k,j,i,1)
 255               x_up(ii)=weitor_up(k,j,i,1)
 256               write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii
 257             enddo
 258             do l=1,nterm_sccor(j,i)
 259               ii=ii+1
 260               ntor_var=ntor_var+1
 261               x(ii)=v2sccor(l,k,j,i)
 262               x_low(ii)=weitor_low(k,j,i,1)
 263               x_up(ii)=weitor_up(k,j,i,1)
 264               write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii
 265             enddo
 266           endif
 267           enddo
 268         enddo
 269       enddo
 270
 271       endif
 272
 273 c Optimized cumulant coefficients
 274
 275       if (mod_fourier(nloctyp).gt.0) then
 276
 277 #ifdef NEWCORR
 278       do i=0,nloctyp-1
 279
 280         if (mod_fourier(i).gt.0) then
 281
 282           do j=1,2
 283             do k=1,3
 284               if (mask_bnew1(k,j,i).gt.0) then
 285                 ii=ii+1
 286                 x(ii)=bnew1(k,j,i)
 287                 x_low(ii)=bnew1_low(k,j,i)
 288                 x_up(ii)=bnew1_up(k,j,i)
 289               endif
 290             enddo
 291           enddo
 292
 293           do j=1,2
 294             do k=1,3
 295               if (mask_bnew2(k,j,i).gt.0) then
 296                 ii=ii+1
 297                 x(ii)=bnew2(k,j,i)
 298                 x_low(ii)=bnew2_low(k,j,i)
 299                 x_up(ii)=bnew2_up(k,j,i)
 300               endif
 301             enddo
 302           enddo
 303
 304           do j=1,2
 305             do k=1,3
 306               if (mask_ccnew(k,j,i).gt.0) then
 307                 ii=ii+1
 308                 x(ii)=ccnew(k,j,i)
 309                 x_low(ii)=ccnew_low(k,j,i)
 310                 x_up(ii)=ccnew_up(k,j,i)
 311               endif
 312             enddo
 313           enddo
 314
 315           do j=1,2
 316             do k=1,3
 317               if (mask_ddnew(k,j,i).gt.0) then
 318                 ii=ii+1
 319                 x(ii)=ddnew(k,j,i)
 320                 x_low(ii)=ddnew_low(k,j,i)
 321                 x_up(ii)=ddnew_up(k,j,i)
 322               endif
 323             enddo
 324           enddo
 325
 326           do j=1,2
 327             if (mask_e0new(j,i).gt.0) then
 328               ii=ii+1
 329               x(ii)=e0new(j,i)
 330               x_low(ii)=e0new_low(j,i)
 331               x_up(ii)=e0new_up(j,i)
 332             endif
 333           enddo
 334
 335           do j=1,2
 336             do k=1,2
 337               do l=1,2
 338                 if (mask_eenew(l,k,j,i).gt.0) then
 339                   ii=ii+1
 340                   x(ii)=eenew(l,k,j,i)
 341                   x_low(ii)=eenew_low(l,k,j,i)
 342                   x_up(ii)=eenew_up(l,k,j,i)
 343                 endif
 344               enddo
 345             enddo
 346           enddo
 347
 348         endif
 349
 350       enddo
 351 #else
 352       do i=1,nloctyp
 353
 354         if (mod_fourier(i).gt.0) then
 355           do j=2,13
 356             if (mask_fourier(j,i).gt.0) then
 357               ii=ii+1
 358               x(ii)=b(j,i)
 359               x_low(ii)=b_low(j,i)
 360               x_up(ii)=b_up(j,i)
 361             endif
 362           enddo
 363         endif
 364
 365       enddo
 366 #endif
 367       endif
 368
 369 c Added SC sigmas and anisotropies
 370
 371       if (mod_side_other) then
 372         do i=1,ntyp
 373           if (mask_sigma(1,i).gt.0) then
 374             ii=ii+1
 375             x(ii)=sigma0(i)
 376             x_low(ii)=sigma_low(1,i)
 377             x_up(ii)=sigma_up(1,i)
 378 c            write (iout,*) "sig0",i,ii,x(ii), x_low(ii),x_up(ii)
 379           endif
 380           if (mask_sigma(2,i).gt.0) then
 381             ii=ii+1
 382             x(ii)=sigii(i)
 383             x_low(ii)=sigma_low(2,i)
 384             x_up(ii)=sigma_up(2,i)
 385 c            write (iout,*) "sigi",i,ii,x(ii), x_low(ii),x_up(ii)
 386           endif
 387           if (mask_sigma(3,i).gt.0) then
 388             ii=ii+1
 389             x(ii)=chip0(i)
 390             x_low(ii)=sigma_low(3,i)
 391             x_up(ii)=sigma_up(3,i)
 392 c            write (iout,*) "chi",i,ii,x(ii), x_low(ii),x_up(ii)
 393           endif
 394           if (mask_sigma(4,i).gt.0) then
 395             ii=ii+1
 396             x(ii)=alp(i)
 397             x_low(ii)=sigma_low(4,i)
 398             x_up(ii)=sigma_up(4,i)
 399 c            write (iout,*) "alp",i,ii,x(ii), x_low(ii),x_up(ii)
 400           endif
 401         enddo
 402       endif
 403
 404
 405       if (mod_elec) then
 406
 407       do i=1,2
 408         do j=1,i
 409           if (mask_elec(i,j,1).gt.0) then
 410             ii=ii+1
 411             x(ii)=epp(i,j)
 412             x_low(ii)=epp_low(i,j)
 413             x_up(ii)=epp_up(i,j)
 414           endif
 415           if (mask_elec(i,j,2).gt.0) then
 416             ii=ii+1
 417             x(ii)=rpp(i,j)
 418             x_low(ii)=rpp_low(i,j)
 419             x_up(ii)=rpp_up(i,j)
 420           endif
 421           if (mask_elec(i,j,3).gt.0) then
 422             ii=ii+1
 423             x(ii)=elpp6(i,j)
 424             x_low(ii)=elpp6_low(i,j)
 425             x_up(ii)=elpp6_up(i,j)
 426           endif
 427           if (mask_elec(i,j,4).gt.0) then
 428             ii=ii+1
 429             x(ii)=elpp3(i,j)
 430             x_low(ii)=elpp3_low(i,j)
 431             x_up(ii)=elpp3_up(i,j)
 432           endif
 433         enddo
 434       enddo
 435
 436       endif
 437 C
 438 C Define the SC-p interaction constants
 439 C
 440       if (mod_scp) then
 441
 442       if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+
 443      &        mask_scp(0,2,2) .gt. 0) then
 444         do j=1,2
 445           if (mask_scp(0,j,1).gt.0) then
 446             ii=ii+1
 447             x(ii)=eps_scp(1,j)
 448             x_low(ii)=epscp_low(0,j)
 449             x_up(ii)=epscp_up(0,j)
 450           endif
 451           if (mask_scp(0,j,2).gt.0) then
 452             ii=ii+1
 453             x(ii)=rscp(1,j)
 454             x_low(ii)=rscp_low(0,j)
 455             x_up(ii)=rscp_up(0,j)
 456           endif
 457         enddo
 458       else
 459         do i=1,20
 460           do j=1,2
 461             if (mask_scp(i,j,1).gt.0) then
 462               ii=ii+1
 463               x(ii)=eps_scp(i,j)
 464               x_low(ii)=epscp_low(i,j)
 465               x_up(ii)=epscp_up(i,j)
 466             endif
 467             if (mask_scp(i,j,2).gt.0) then
 468               ii=ii+1
 469               x(ii)=rscp(i,j)
 470               x_low(ii)=rscp_low(i,j)
 471               x_up(ii)=rscp_up(i,j)
 472             endif
 473           enddo
 474         enddo
 475       endif
 476
 477       endif
 478
 479
 480 c Compute boundaries, if defined relative to starting values
 481       do i=1,ii
 482         write (iout,*) "i",i," x_low",x_low(i)," x_up",x_up(i)
 483         if (x_up(i).eq.0.0d0 .and. x_low(i).gt.0.0d0) then
 484           xchuj=x_low(i)
 485           x_low(i)=x(i)-xchuj
 486           x_up(i)=x(i)+xchuj
 487         else if (x_up(i).eq.-1.0d0 .and. x_low(i).gt.0.0d0) then
 488           xchuj=x_low(i)
 489           if (x(i).gt.0.0d0) then
 490             x_low(i)=x(i)*(1.0d0-xchuj)
 491             x_up(i)=x(i)*(1.0d0+xchuj)
 492           else
 493             x_low(i)=x(i)*(1.0d0+xchuj)
 494             x_up(i)=x(i)*(1.0d0-xchuj)
 495           endif
 496         endif
 497       enddo
 498
 499 c Base of heat capacity
 500
 501       do iprot=1,nprot
 502         if (nbeta(2,iprot).gt.0) then
 503           ii=ii+1
 504           x(ii)=heatbase(iprot)
 505           x_low(ii)=-10.0d0
 506           x_up(ii)=10.0d0
 507         endif
 508       enddo
 509
 510       nvarr=ii
 511
 512       write (iout,*) "Number of optimizable parameters:",nvarr
 513       write (iout,*)
 514       write (iout,*) "Initial variables and boundaries:"
 515       write (iout,'(3x,3a10)') "x","x_low","x_up"
 516       do i=1,nvarr
 517         write (iout,'(i3,3f10.5)') i,x(i),x_low(i),x_up(i)
 518       enddo
 519       call flush(iout)
 520       return
 521       end
 522 c------------------------------------------------------------------------------
 523       subroutine x2w(nvarr,x)
 524       implicit none
 525       include "DIMENSIONS"
 526       include "DIMENSIONS.ZSCOPT"
 527       include "COMMON.CONTROL"
 528       include "COMMON.NAMES"
 529       include "COMMON.WEIGHTS"
 530       include "COMMON.INTERACT"
 531       include "COMMON.ENERGIES"
 532       include "COMMON.FFIELD"
 533       include "COMMON.TORSION"
 534       include "COMMON.IOUNITS"
 535       include "COMMON.VMCPAR"
 536       include "COMMON.CLASSES"
 537       include "COMMON.WEIGHTDER"
 538       include "COMMON.SCCOR"
 539       logical lprint
 540       integer i,j,k,kk,l,ll,ii,nvarr,iti,itj,iprot,iblock
 541       double precision rri,v0ij,si
 542       double precision x(nvarr)
 543       double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2,
 544      &  ratsig1,ratsig2,rsum_max,dwa16,r_augm
 545       double precision eps_temp(ntyp,ntyp)
 546       double precision ftune_eps,ftune_epsprim
 547       logical in_sc_pair_list
 548       external in_sc_pair_list
 549       integer itor2typ(3) /10,9,20/
 550
 551       lprint=.false.
 552
 553       if (lprint) write (iout,*) "x2w: nvar",nvarr
 554       if (lprint) write (iout,*) "x",(x(i),i=1,nvarr)
 555       ii=0
 556       do i=1,n_ene
 557         if (imask(i).eq.1) then
 558           ii=ii+1
 559           ww(i)=x(ii)
 560         endif
 561       enddo
 562
 563       wsc=ww(1)
 564       wscp=ww(2)
 565       welec=ww(3)
 566       wcorr=ww(4)
 567       wcorr5=ww(5)
 568       wcorr6=ww(6)
 569       wel_loc=ww(7)
 570       wturn3=ww(8)
 571       wturn4=ww(9)
 572       wturn6=ww(10)
 573       wang=ww(11)
 574       wscloc=ww(12)
 575       wtor=ww(13)
 576       wtor_d=ww(14)
 577       wstrain=ww(15)
 578       wvdwpp=ww(16)
 579       wbond=ww(17)
 580       wsccor=ww(19)
 581 #ifdef SCP14
 582       scal14=ww(18)/ww(2)
 583 #endif
 584       do i=1,n_ene
 585         weights(i)=ww(i)
 586       enddo
 587
 588       if (lprint) then
 589         write (iout,*) "Weights:"
 590         do i=1,n_ene
 591           write (iout,'(a,f10.5)') wname(i),ww(i)
 592         enddo
 593       endif
 594
 595 C SC parameters
 596       if (mod_side) then
 597         do i=1,nsingle_sc
 598           ii=ii+1
 599           iti=ityp_ssc(i)
 600           eps(iti,iti)=x(ii)
 601           do j=1,iti-1
 602             if (.not. in_sc_pair_list(iti,j)) then
 603               eps(iti,j)=eps_orig(iti,j)
 604      &           +0.5d0*(x(ii)-eps_orig(iti,iti))
 605               eps(j,iti)=eps(iti,j)
 606             endif
 607           enddo
 608         enddo
 609         do i=1,npair_sc
 610           ii=ii+1
 611           iti=ityp_psc(1,i)
 612           itj=ityp_psc(2,i)
 613           eps(iti,itj)=x(ii)
 614           eps(itj,iti)=x(ii)
 615         enddo
 616       endif
 617
 618 #ifdef NEWCORR
 619       if (mod_tor .or. mod_fourier(nloctyp).gt.0) then
 620 c      write (iout,*) "w2x: mod_tor",mod_tor," mod_fourier",
 621 c     &   mod_fourier(nloctyp)
 622 #else
 623       if (mod_tor) then
 624 #endif
 625       do i=-ntortyp+1,ntortyp-1
 626         do j=-ntortyp+1,ntortyp-1
 627           if (tor_mode.eq.0) then
 628           do iblock=1,2
 629             if (mask_tor(0,j,i,iblock).eq.1) then
 630               ii=ii+1
 631               weitor(0,j,i,iblock)=x(ii)
 632             else if (mask_tor(0,j,i,iblock).eq.2) then
 633               do k=1,nterm(j,i,iblock)
 634                 ii=ii+1
 635                 v1(k,j,i,iblock)=x(ii)
 636               enddo
 637             else if (mask_tor(0,j,i,iblock).gt.2) then
 638               do k=1,nterm(j,i,iblock)
 639                 ii=ii+1
 640                 v1(k,j,i,iblock)=x(ii)
 641               enddo
 642               do k=1,nterm(j,i,iblock)
 643                 ii=ii+1
 644                 v2(k,j,i,iblock)=x(ii)
 645               enddo
 646             endif
 647             v0ij=0.0d0
 648             si=-1.0d0
 649             do k=1,nterm(i,j,iblock)
 650               v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 651               v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 652               v0ij=v0ij+si*v1(k,i,j,iblock)
 653               si=-si
 654             enddo
 655             do k=1,nlor(i,j,iblock)
 656               v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 657             enddo
 658             v0(i,j,iblock)=v0ij
 659             v0(-i,-j,iblock)=v0ij
 660           enddo
 661           else if (tor_mode.eq.1) then
 662             if (mask_tor(0,j,i,1).eq.1) then
 663               ii=ii+1
 664               weitor(0,j,i,1)=x(ii)
 665             else if (mask_tor(0,j,i,1).eq.2) then
 666               do k=1,nterm_kcc(j,i)
 667                 do l=1,nterm_kcc_Tb(j,i)
 668                   do ll=1,nterm_kcc_Tb(j,i)
 669                     ii=ii+1
 670                     v1_kcc(ll,l,k,j,i)=x(ii)
 671                   enddo ! ll
 672                 enddo ! l
 673               enddo
 674             else if (mask_tor(0,j,i,1).gt.2) then
 675               do k=1,nterm_kcc(j,i)
 676                 do l=1,nterm_kcc_Tb(j,i)
 677                   do ll=1,nterm_kcc_Tb(j,i)
 678                     ii=ii+1
 679                     v1_kcc(ll,l,k,j,i)=x(ii)
 680                   enddo ! ll
 681                 enddo ! l
 682               enddo
 683               do k=1,nterm_kcc(j,i)
 684                 do l=1,nterm_kcc_Tb(j,i)
 685                   do ll=1,nterm_kcc_Tb(j,i)
 686                     ii=ii+1
 687                     v2_kcc(ll,l,k,j,i)=x(ii)
 688                   enddo! ll
 689                 enddo ! l
 690               enddo
 691             endif
 692           else
 693 #ifdef NEWCORR
 694 c Compute torsional coefficients from local energy surface expansion coeffs
 695             if (mask_tor(0,j,i,1).eq.1) then
 696               ii=ii+1
 697               weitor(0,j,i,1)=x(ii)
 698 c            else if (mask_tor(0,j,i,1).eq.2) then
 699             endif
 700 #else
 701           else
 702             write (iout,*) "ERROR: TOR_MODE>1 ALLOWED ONLY WITH NEWCORR"
 703             stop
 704 #endif
 705           endif
 706         enddo
 707       enddo
 708
 709       endif
 710       if (mod_sccor) then
 711
 712       do i=-nsccortyp,nsccortyp
 713         do j=-nsccortyp,nsccortyp
 714           do k=1,3
 715           if (mask_tor(k,j,i,1).eq.1) then
 716             ii=ii+1
 717             weitor(k,j,i,1)=x(ii)
 718           else if (mask_tor(k,j,i,1).eq.2) then
 719             do l=1,nterm(j,i,1)
 720               ii=ii+1
 721               v1sccor(l,k,j,i)=x(ii)
 722             enddo
 723           else if (mask_tor(k,j,i,1).gt.2) then
 724             do l=1,nterm_sccor(j,i)
 725               ii=ii+1
 726               v1sccor(l,k,j,i)=x(ii)
 727             enddo
 728             do l=1,nterm_sccor(j,i)
 729               ii=ii+1
 730               v2sccor(l,k,j,i)=x(ii)
 731             enddo
 732           endif
 733           enddo
 734         enddo
 735       enddo
 736
 737       endif
 738
 739       if (mod_fourier(nloctyp).gt.0) then
 740
 741 #ifdef NEWCORR
 742       do i=0,nloctyp-1
 743
 744         if (mod_fourier(i).gt.0) then
 745
 746           do j=1,2
 747             do k=1,3
 748               if (mask_bnew1(k,j,i).gt.0) then
 749                 ii=ii+1
 750                 bnew1(k,j,i)=x(ii)
 751               endif
 752             enddo
 753           enddo
 754
 755           do j=1,2
 756             do k=1,3
 757               if (mask_bnew2(k,j,i).gt.0) then
 758                 ii=ii+1
 759                 bnew2(k,j,i)=x(ii)
 760               endif
 761             enddo
 762           enddo
 763
 764           do j=1,2
 765             do k=1,3
 766               if (mask_ccnew(k,j,i).gt.0) then
 767                 ii=ii+1
 768                 ccnew(k,j,i)=x(ii)
 769               endif
 770             enddo
 771           enddo
 772
 773           do j=1,2
 774             do k=1,3
 775               if (mask_ddnew(k,j,i).gt.0) then
 776                 ii=ii+1
 777                 ddnew(k,j,i)=x(ii)
 778               endif
 779             enddo
 780           enddo
 781
 782           do j=1,2
 783             if (mask_e0new(j,i).gt.0) then
 784               ii=ii+1
 785               e0new(j,i)=x(ii)
 786             endif
 787           enddo
 788
 789           do j=1,2
 790             do k=1,2
 791               do l=1,2
 792                 if (mask_eenew(l,k,j,i).gt.0) then
 793                   ii=ii+1
 794                   eenew(l,k,j,i)=x(ii)
 795                 endif
 796               enddo
 797             enddo
 798           enddo
 799
 800         endif
 801
 802       enddo
 803
 804       do i=1,nloctyp-1
 805         do j=1,3
 806           bnew1(j,1,-i)= bnew1(j,1,i)
 807           bnew1(j,2,-i)=-bnew1(j,2,i)
 808           bnew2(j,1,-i)= bnew2(j,1,i)
 809           bnew2(j,2,-i)=-bnew2(j,2,i)
 810         enddo
 811         do j=1,3
 812 c          ccnew(j,1,i)=ccnew(j,1,i)/2
 813 c          ccnew(j,2,i)=ccnew(j,2,i)/2
 814 c          ddnew(j,1,i)=ddnew(j,1,i)/2
 815 c          ddnew(j,2,i)=ddnew(j,2,i)/2
 816           ccnew(j,1,-i)=ccnew(j,1,i)
 817           ccnew(j,2,-i)=-ccnew(j,2,i)
 818           ddnew(j,1,-i)=ddnew(j,1,i)
 819           ddnew(j,2,-i)=-ddnew(j,2,i)
 820         enddo
 821         e0new(1,-i)= e0new(1,i)
 822         e0new(2,-i)=-e0new(2,i)
 823         e0new(3,-i)=-e0new(3,i)
 824         do kk=1,2
 825           eenew(kk,1,1,-i)= eenew(kk,1,1,i)
 826           eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
 827           eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
 828           eenew(kk,2,2,-i)= eenew(kk,2,2,i)
 829         enddo
 830       enddo
 831       if (lprint) then
 832         write (iout,'(a)') "Coefficients of the multibody terms"
 833         do i=-nloctyp+1,nloctyp-1
 834           write (iout,*) "Type: ",onelet(iloctyp(i))
 835           write (iout,*) "Coefficients of the expansion of B1"
 836           do j=1,2
 837             write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
 838           enddo
 839           write (iout,*) "Coefficients of the expansion of B2"
 840           do j=1,2
 841             write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
 842           enddo
 843           write (iout,*) "Coefficients of the expansion of C"
 844           write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
 845           write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
 846           write (iout,*) "Coefficients of the expansion of D"
 847           write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
 848           write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
 849           write (iout,*) "Coefficients of the expansion of E"
 850           write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
 851           do j=1,2
 852             do k=1,2
 853               write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
 854             enddo
 855           enddo
 856         enddo
 857         call flush(iout)
 858       endif
 859       if (tor_mode.eq.2) then
 860 c Recalculate new torsional potential coefficients
 861         do i=-ntortyp+1,ntortyp-1
 862           do j=-ntortyp+1,ntortyp-1
 863             do k=1,nterm_kcc_Tb(j,i)
 864               do l=1,nterm_kcc_Tb(j,i)
 865                 v1_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,1,j)
 866      &                         +bnew1(k,2,i)*bnew2(l,2,j)
 867                 v2_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,2,j)
 868      &                         +bnew1(k,2,i)*bnew2(l,1,j)
 869               enddo
 870             enddo
 871             do k=1,nterm_kcc_Tb(j,i)
 872               do l=1,nterm_kcc_Tb(j,i)
 873                 v1_kcc(k,l,2,i,j)=-0.25*(ccnew(k,1,i)*ddnew(l,1,j)
 874      &                         -ccnew(k,2,i)*ddnew(l,2,j))
 875                 v2_kcc(k,l,2,i,j)=-0.25*(ccnew(k,2,i)*ddnew(l,1,j)
 876      &                         +ccnew(k,1,i)*ddnew(l,2,j))
 877               enddo
 878             enddo
 879           enddo
 880         enddo
 881       endif
 882 #else
 883       do i=0,nloctyp-1
 884
 885         if (mod_fourier(i).gt.0) then
 886
 887           do j=2,13
 888             if (mask_fourier(j,i).gt.0) then
 889               ii=ii+1
 890               b(j,i)=x(ii)
 891             endif
 892           enddo
 893
 894           if (lprint) then
 895             write (iout,*) "Fourier coefficients of cumulants"
 896             write (iout,*) 'Type',i
 897             write (iout,'(a,i2,a,f10.5)') ('b(',j,')=',b(j,i),j=1,13)
 898           endif
 899
 900           if (i.gt.0) then
 901             b(2,-i)= b(2,i)
 902             b(3,-i)= b(3,i)
 903             b(4,-i)=-b(4,i)
 904             b(5,-i)=-b(5,i)
 905           endif
 906 c          B1(1,i)  = b(3,i)
 907 c          B1(2,i)  = b(5,i)
 908 c          B1tilde(1,i) = b(3,i)
 909 c          B1tilde(2,i) =-b(5,i)
 910 c          B2(1,i)  = b(2,i)
 911 c          B2(2,i)  = b(4,i)
 912           CCold(1,1,i)= b(7,i)
 913           CCold(2,2,i)=-b(7,i)
 914           CCold(2,1,i)= b(9,i)
 915           CCold(1,2,i)= b(9,i)
 916           CCold(1,1,-i)= b(7,i)
 917           CCold(2,2,-i)=-b(7,i)
 918           CCold(2,1,-i)=-b(9,i)
 919           CCold(1,2,-i)=-b(9,i)
 920 c          Ctilde(1,1,i)=b(7,i)
 921 c          Ctilde(1,2,i)=b(9,i)
 922 c          Ctilde(2,1,i)=-b(9,i)
 923 c          Ctilde(2,2,i)=b(7,i)
 924           DDold(1,1,i)= b(6,i)
 925           DDold(2,2,i)=-b(6,i)
 926           DDold(2,1,i)= b(8,i)
 927           DDold(1,2,i)= b(8,i)
 928           DDold(1,1,-i)= b(6,i)
 929           DDold(2,2,-i)=-b(6,i)
 930           DDold(2,1,-i)=-b(8,i)
 931           DDold(1,2,-i)=-b(8,i)
 932 c          Dtilde(1,1,i)=b(6,i)
 933 c          Dtilde(1,2,i)=b(8,i)
 934 c          Dtilde(2,1,i)=-b(8,i)
 935 c          Dtilde(2,2,i)=b(6,i)
 936           EEold(1,1,i)= b(10,i)+b(11,i)
 937           EEold(2,2,i)=-b(10,i)+b(11,i)
 938           EEold(2,1,i)= b(12,i)-b(13,i)
 939           EEold(1,2,i)= b(12,i)+b(13,i)
 940           EEold(1,1,i)= b(10,i)+b(11,i)
 941           EEold(2,2,i)=-b(10,i)+b(11,i)
 942           EEold(2,1,i)= b(12,i)-b(13,i)
 943           EEold(1,2,i)= b(12,i)+b(13,i)
 944           EEold(1,1,-i)= b(10,i)+b(11,i)
 945           EEold(2,2,-i)=-b(10,i)+b(11,i)
 946           EEold(2,1,-i)=-b(12,i)+b(13,i)
 947           EEold(1,2,-i)=-b(12,i)-b(13,i)
 948
 949         endif
 950
 951       enddo
 952       if (lprint) then
 953       write (iout,*)
 954       write (iout,*)
 955      &"Coefficients of the cumulants (independent of valence angles)"
 956       do i=-nloctyp+1,nloctyp-1
 957         write (iout,*) 'Type ',onelet(iloctyp(i))
 958         write (iout,*) 'B1'
 959         write(iout,'(2f10.5)') B(3,i),B(5,i)
 960         write (iout,*) 'B2'
 961         write(iout,'(2f10.5)') B(2,i),B(4,i)
 962         write (iout,*) 'CC'
 963         do j=1,2
 964           write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
 965         enddo
 966         write(iout,*) 'DD'
 967         do j=1,2
 968           write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
 969         enddo
 970         write(iout,*) 'EE'
 971         do j=1,2
 972           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
 973         enddo
 974       enddo
 975       endif
 976 #endif
 977
 978       endif
 979 c SC sigmas and anisotropies as parameters
 980       if (mod_side_other) then
 981 c        do i=1,ntyp
 982 c          write (iout,*) "mask",i,(mask_sigma(k,i),k=1,4)
 983 c        enddo
 984         do i=1,ntyp
 985           if (mask_sigma(1,i).gt.0) then
 986             ii=ii+1
 987             sigma0(i)=x(ii)
 988 c            write(iout,*) "sig0",i,ii,x(ii)
 989           endif
 990           if (mask_sigma(2,i).gt.0) then
 991             ii=ii+1
 992             sigii(i)=x(ii)
 993 c            write(iout,*) "sigi",i,ii,x(ii)
 994           endif
 995           if (mask_sigma(3,i).gt.0) then
 996             ii=ii+1
 997             chip0(i)=x(ii)
 998 c            write(iout,*) "chip",i,ii,x(ii)
 999           endif
1000           if (mask_sigma(4,i).gt.0) then
1001             ii=ii+1
1002             alp(i)=x(ii)
1003 c            write(iout,*) "alp",i,ii,x(ii)
1004           endif
1005         enddo
1006       endif
1007
1008       if (mod_side.or.mod_side_other) then
1009         if (lprint) write(iout,'(/3a,2i3)') 'Potential is ',
1010      &    potname(ipot),', exponents are ',expon,2*expon
1011 C-----------------------------------------------------------------------
1012 C Calculate the "working" parameters of SC interactions.
1013         dwa16=2.0d0**(1.0d0/6.0d0)
1014         if (ipot.eq.4) then
1015           do i=1,ntyp
1016             chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1017           enddo
1018         endif
1019         do i=1,ntyp
1020           do j=i,ntyp
1021             sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1022             sigma(j,i)=sigma(i,j)
1023             rs0(i,j)=dwa16*sigma(i,j)
1024             rs0(j,i)=rs0(i,j)
1025           enddo
1026         enddo
1027         if (lprint) call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1028         if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1029      &   'Working parameters of the SC interactions:',
1030      &   '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1031      &   '  chi1   ','   chi2   '
1032         do i=1,ntyp
1033           do j=i,ntyp
1034             epsij=eps(i,j)
1035             if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1036               rrij=sigma(i,j)
1037             else
1038               rrij=rr0(i)+rr0(j)
1039             endif
1040             r0(i,j)=rrij
1041             r0(j,i)=rrij
1042             rrij=rrij**expon
1043             epsij=eps(i,j)
1044 c           sigeps=dsign(1.0D0,epsij)
1045             aa(i,j)=ftune_eps(epsij)*rrij*rrij
1046             bb(i,j)=-epsij*rrij
1047             aa(j,i)=aa(i,j)
1048             bb(j,i)=bb(i,j)
1049             if (ipot.gt.2) then
1050               sigt1sq=sigma0(i)**2
1051               sigt2sq=sigma0(j)**2
1052               sigii1=sigii(i)
1053               sigii2=sigii(j)
1054               ratsig1=sigt2sq/sigt1sq
1055               ratsig2=1.0D0/ratsig1
1056               chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1057               if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1058               rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1059             else
1060               rsum_max=sigma(i,j)
1061             endif
1062             sigmaii(i,j)=rsum_max
1063             sigmaii(j,i)=rsum_max
1064             if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max)
1065      &      then
1066               r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1067               augm(i,j)=epsij*r_augm**(2*expon)
1068 c             augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1069               augm(j,i)=augm(i,j)
1070             else
1071               augm(i,j)=0.0D0
1072               augm(j,i)=0.0D0
1073             endif
1074             if (lprint) then
1075               write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1076      &        restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1077      &        sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1078             endif
1079           enddo
1080         enddo
1081       endif
1082
1083
1084       if (mod_elec) then
1085
1086       do i=1,2
1087         do j=1,i
1088           if (mask_elec(i,j,1).gt.0) then
1089             ii=ii+1
1090             epp(i,j)=x(ii)
1091             if (epp(i,j).lt.0.0d0) epp(i,j)=0.0d0
1092             epp(j,i)=x(ii)
1093           endif
1094           if (mask_elec(i,j,2).gt.0) then
1095             ii=ii+1
1096             rpp(i,j)=dabs(x(ii))
1097             rpp(j,i)=dabs(x(ii))
1098           endif
1099           if (mask_elec(i,j,3).gt.0) then
1100             ii=ii+1
1101             elpp6(i,j)=x(ii)
1102             if (elpp6(i,j).gt.0.0d0) elpp6(i,j)=0.0d0
1103             elpp6(j,i)=elpp6(i,j)
1104           endif
1105           if (mask_elec(i,j,4).gt.0) then
1106             ii=ii+1
1107             elpp3(i,j)=x(ii)
1108             elpp3(j,i)=x(ii)
1109           endif
1110           rri=rpp(i,j)**6
1111           app (i,j)=epp(i,j)*rri*rri
1112           app (j,i)=epp(j,i)*rri*rri
1113           bpp (i,j)=-2.0D0*epp(i,j)*rri
1114           bpp (j,i)=-2.0D0*epp(i,j)*rri
1115           ael6(i,j)=elpp6(i,j)*4.2D0**6
1116           ael6(j,i)=elpp6(i,j)*4.2D0**6
1117           ael3(i,j)=elpp3(i,j)*4.2D0**3
1118           ael3(j,i)=elpp3(i,j)*4.2D0**3
1119         enddo
1120       enddo
1121       if (lprint) then
1122         write (iout,'(/a)') 'Electrostatic interaction constants:'
1123         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1124      &            'IT','JT','APP','BPP','AEL6','AEL3'
1125         do i=1,2
1126           do j=1,2
1127             write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1128      &                    ael6(i,j),ael3(i,j)
1129           enddo
1130         enddo
1131       endif
1132
1133       endif
1134 C
1135 C Define the SC-p interaction constants
1136 C
1137       if (mod_scp) then
1138
1139       if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+
1140      &        mask_scp(0,2,2) .gt. 0) then
1141         do j=1,2
1142           if (mask_scp(0,j,1).gt.0) then
1143             ii=ii+1
1144             do i=1,20
1145               eps_scp(i,j)=x(ii)
1146             enddo
1147           endif
1148           if (mask_scp(0,j,2).gt.0) then
1149             ii=ii+1
1150             do i=1,20
1151               rscp(i,j)=x(ii)
1152             enddo
1153           endif
1154         enddo
1155       else
1156         do i=1,20
1157           do j=1,2
1158             if (mask_scp(i,j,1).gt.0) then
1159               ii=ii+1
1160               eps_scp(i,j)=x(ii)
1161             endif
1162             if (mask_scp(i,j,2).gt.0) then
1163               ii=ii+1
1164               rscp(i,j)=x(ii)
1165             endif
1166           enddo
1167         enddo
1168       endif
1169       do i=1,20
1170         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1171         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1172         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1173         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1174       enddo
1175
1176       if (lprint) then
1177         write (iout,*) "Parameters of SC-p interactions:"
1178         do i=1,20
1179           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1180      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1181         enddo
1182       endif
1183
1184       endif
1185
1186
1187 c      write (iout,*) "x2w: nvar",nvarr
1188
1189 c Base of heat capacity
1190
1191       do iprot=1,nprot
1192         if (nbeta(2,iprot).gt.0) then
1193           ii=ii+1
1194           heatbase(iprot)=x(ii)
1195         endif
1196       enddo
1197
1198       if (ii.ne.nvarr) then
1199         write (iout,*) "!!!!! SEVERE ERROR IN X2W !!!!! (chuj nastapil)"
1200         write (iout,*) "NVAR=",nvarr," LAST VARIABLE INDEX=",ii
1201         stop "CHUJ ABSOLUTNY!!!!"
1202       endif
1203
1204       return
1205       end
1206 c---------------------------------------------------------------------------
1207       double precision function ftune_eps(x)
1208       implicit none
1209       double precision x,y
1210       double precision delta /1.0d-5/
1211       if (dabs(x).lt.delta) then
1212         y=x/delta
1213         ftune_eps=delta*(3.0d0+6.0d0*y**2-y**4)/8
1214       else
1215         ftune_eps=dabs(x)
1216       endif
1217       return
1218       end
1219 c---------------------------------------------------------------------------
1220       double precision function ftune_epsprim(x)
1221       implicit none
1222       double precision x,y
1223       double precision delta /1.0d-5/
1224       if (dabs(x).lt.delta) then
1225         y=x/delta
1226         ftune_epsprim=1.5d0*y-0.5*y**3
1227       else
1228         ftune_epsprim=dsign(1.0d0,x)
1229       endif
1230       return
1231       end
1232 c---------------------------------------------------------------------------
1233       double precision function ftune_epsbis(x)
1234       implicit none
1235       double precision x,y
1236       double precision delta /1.0d-5/
1237       if (dabs(x).lt.delta) then
1238         y=x/delta
1239         ftune_epsbis=(1.5d0-1.5d0*y**3)/delta
1240       else
1241         ftune_epsbis=0.0d0
1242       endif
1243       return
1244       end
1245 c---------------------------------------------------------------------------
1246       logical function in_sc_pair_list(iti,itj)
1247       implicit none
1248       include "DIMENSIONS"
1249       include "DIMENSIONS.ZSCOPT"
1250       include "COMMON.WEIGHTS"
1251       integer i,iti,itj
1252       do i=1,npair_sc
1253         if (ityp_psc(1,i).eq.iti .and. ityp_psc(2,i).eq.itj .or.
1254      &      ityp_psc(1,i).eq.itj .and. ityp_psc(2,i).eq.iti) then
1255           in_sc_pair_list=.true.
1256           return
1257         endif
1258       enddo
1259       in_sc_pair_list=.false.
1260       return
1261       end