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