source/maxlik/src_FPy.org/w2x_eps.F

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