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