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

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