source/maxlik/src_MD_T_maxlik-NEWCORR-PMF/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         enddo
 608       endif
 609
 610
 611       if (mod_elec) then
 612
 613       do i=1,2
 614         do j=1,i
 615           if (mask_elec(i,j,1).gt.0) then
 616             ii=ii+1
 617             x(ii)=epp(i,j)
 618             x_low(ii)=epp_low(i,j)
 619             x_up(ii)=epp_up(i,j)
 620           endif
 621           if (mask_elec(i,j,2).gt.0) then
 622             ii=ii+1
 623             x(ii)=rpp(i,j)
 624             x_low(ii)=rpp_low(i,j)
 625             x_up(ii)=rpp_up(i,j)
 626           endif
 627           if (mask_elec(i,j,3).gt.0) then
 628             ii=ii+1
 629             x(ii)=elpp6(i,j)
 630             x_low(ii)=elpp6_low(i,j)
 631             x_up(ii)=elpp6_up(i,j)
 632           endif
 633           if (mask_elec(i,j,4).gt.0) then
 634             ii=ii+1
 635             x(ii)=elpp3(i,j)
 636             x_low(ii)=elpp3_low(i,j)
 637             x_up(ii)=elpp3_up(i,j)
 638           endif
 639         enddo
 640       enddo
 641
 642       endif
 643 C
 644 C Define the SC-p interaction constants
 645 C
 646       if (mod_scp) then
 647
 648       if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+
 649      &        mask_scp(0,2,2) .gt. 0) then
 650         do j=1,2
 651           if (mask_scp(0,j,1).gt.0) then
 652             ii=ii+1
 653             x(ii)=eps_scp(1,j)
 654             x_low(ii)=epscp_low(0,j)
 655             x_up(ii)=epscp_up(0,j)
 656           endif
 657           if (mask_scp(0,j,2).gt.0) then
 658             ii=ii+1
 659             x(ii)=rscp(1,j)
 660             x_low(ii)=rscp_low(0,j)
 661             x_up(ii)=rscp_up(0,j)
 662           endif
 663         enddo
 664       else
 665         do i=1,20
 666           do j=1,2
 667             if (mask_scp(i,j,1).gt.0) then
 668               ii=ii+1
 669               x(ii)=eps_scp(i,j)
 670               x_low(ii)=epscp_low(i,j)
 671               x_up(ii)=epscp_up(i,j)
 672             endif
 673             if (mask_scp(i,j,2).gt.0) then
 674               ii=ii+1
 675               x(ii)=rscp(i,j)
 676               x_low(ii)=rscp_low(i,j)
 677               x_up(ii)=rscp_up(i,j)
 678             endif
 679           enddo
 680         enddo
 681       endif
 682
 683       endif
 684
 685 c Weight of the shielding term
 686       if (imask(ind_shield).eq.1) then
 687         ii=ii+1
 688         x(ii)=ww(ind_shield)
 689         x_low(ii)=ww_low(ind_shield)
 690         x_up(ii)=ww_up(ind_shield)
 691       endif
 692
 693
 694 c Compute boundaries, if defined relative to starting values
 695       do i=1,ii
 696         write (iout,*) "i",i," x_low",x_low(i)," x_up",x_up(i)
 697         if (x_up(i).eq.0.0d0 .and. x_low(i).gt.0.0d0) then
 698           xchuj=x_low(i)
 699           x_low(i)=x(i)-xchuj
 700           x_up(i)=x(i)+xchuj
 701         else if (x_up(i).eq.-1.0d0 .and. x_low(i).gt.0.0d0) then
 702           xchuj=x_low(i)
 703           if (x(i).gt.0.0d0) then
 704             x_low(i)=x(i)*(1.0d0-xchuj)
 705             x_up(i)=x(i)*(1.0d0+xchuj)
 706           else
 707             x_low(i)=x(i)*(1.0d0+xchuj)
 708             x_up(i)=x(i)*(1.0d0-xchuj)
 709           endif
 710         endif
 711       enddo
 712
 713 c Base of heat capacity
 714
 715       do iprot=1,nprot
 716         if (nbeta(2,iprot).gt.0) then
 717           ii=ii+1
 718           x(ii)=heatbase(iprot)
 719           x_low(ii)=-10.0d0
 720           x_up(ii)=10.0d0
 721         endif
 722       enddo
 723
 724       nvarr=ii
 725
 726       write (iout,*) "Number of optimizable parameters:",nvarr
 727       write (iout,*)
 728       write (iout,*) "Initial variables and boundaries:"
 729       write (iout,'(3x,3a15)') "     x","     x_low","     x_up"
 730       do i=1,nvarr
 731         write (iout,'(i3,3e15.5)') i,x(i),x_low(i),x_up(i)
 732       enddo
 733       call flush(iout)
 734       return
 735       end
 736 c------------------------------------------------------------------------------
 737       subroutine x2w(nvarr,x)
 738       implicit none
 739       include "DIMENSIONS"
 740       include "DIMENSIONS.ZSCOPT"
 741       include "COMMON.CONTROL"
 742       include "COMMON.NAMES"
 743       include "COMMON.WEIGHTS"
 744       include "COMMON.INTERACT"
 745       include "COMMON.ENERGIES"
 746       include "COMMON.FFIELD"
 747       include "COMMON.TORSION"
 748       include "COMMON.LOCAL"
 749       include "COMMON.IOUNITS"
 750       include "COMMON.VMCPAR"
 751       include "COMMON.CLASSES"
 752       include "COMMON.WEIGHTDER"
 753       include "COMMON.SCCOR"
 754       include "COMMON.SHIELD"
 755       logical lprint
 756       integer i,j,k,kk,l,ll,ii,nvarr,iti,itj,iprot,iblock
 757       double precision rri,v0ij,si
 758       double precision x(nvarr)
 759       double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2,
 760      &  ratsig1,ratsig2,rsum_max,dwa16,r_augm
 761       double precision eps_temp(ntyp,ntyp)
 762       double precision ftune_eps,ftune_epsprim
 763       logical in_sc_pair_list
 764       external in_sc_pair_list
 765       integer itor2typ(3) /10,9,20/
 766       integer ind_shield /25/
 767
 768       lprint=.false.
 769
 770       if (lprint) write (iout,*) "x2w: nvar",nvarr
 771       if (lprint) write (iout,*) "x",(x(i),i=1,nvarr)
 772       if (lprint) write(iout,*)"X2W: SPLIT_FOURIERTOR",SPLIT_FOURIERTOR,
 773      & " tor_mode",tor_mode
 774       ii=0
 775       do i=1,n_ene
 776         if (imask(i).eq.1 .and. i.ne.ind_shield) then
 777           ii=ii+1
 778           ww(i)=x(ii)
 779 c          write (iout,*) i, wname(i)," ",ww(i)
 780         endif
 781       enddo
 782 #ifdef NEWCORR
 783 c 2/10/208 AL: insert the base torsional PMFs
 784       if (mod_fourier(nloctyp).gt.0) then
 785         do iti=0,2
 786           do itj=-2,2
 787             if (mask_e0(iti,itj).gt.0) then
 788               ii=ii+1
 789               e0(iti,itj)=x(ii)
 790             endif
 791           enddo
 792         enddo
 793         if (mask_ello_PMF.gt.0) then
 794           ii=ii+1
 795           wello_PMF=x(ii)
 796         endif
 797         if (mask_turn_PMF.gt.0) then
 798           ii=ii+1
 799           wturn_PMF=x(ii)
 800         endif
 801       endif
 802 #endif
 803       wsc=ww(1)
 804       wscp=ww(2)
 805       welec=ww(3)
 806       wcorr=ww(4)
 807       wcorr5=ww(5)
 808       wcorr6=ww(6)
 809       wel_loc=ww(7)
 810       wturn3=ww(8)
 811       wturn4=ww(9)
 812       wturn6=ww(10)
 813       wang=ww(11)
 814       wscloc=ww(12)
 815       wtor=ww(13)
 816       wtor_d=ww(14)
 817       wstrain=ww(15)
 818       wvdwpp=ww(16)
 819       wbond=ww(17)
 820       wsccor=ww(19)
 821 #ifdef SCP14
 822       scal14=ww(18)/ww(2)
 823 #endif
 824       do i=1,n_ene
 825         weights(i)=ww(i)
 826       enddo
 827
 828       if (lprint) then
 829         write (iout,*) "Weights:"
 830         do i=1,n_ene
 831           write (iout,'(a,f10.5)') wname(i),ww(i)
 832         enddo
 833       endif
 834
 835 C SC parameters
 836       if (mod_side) then
 837         do i=1,nsingle_sc
 838           ii=ii+1
 839           iti=ityp_ssc(i)
 840           eps(iti,iti)=x(ii)
 841           do j=1,iti-1
 842             if (.not. in_sc_pair_list(iti,j)) then
 843               eps(iti,j)=eps_orig(iti,j)
 844      &           +0.5d0*(x(ii)-eps_orig(iti,iti))
 845               eps(j,iti)=eps(iti,j)
 846             endif
 847           enddo
 848         enddo
 849         do i=1,npair_sc
 850           ii=ii+1
 851           iti=ityp_psc(1,i)
 852           itj=ityp_psc(2,i)
 853           eps(iti,itj)=x(ii)
 854           eps(itj,iti)=x(ii)
 855         enddo
 856       endif
 857
 858 #ifdef NEWCORR
 859       if (mod_tor .or. mod_fourier(nloctyp).gt.0) then
 860 c      write (iout,*) "w2x: mod_tor",mod_tor," mod_fourier",
 861 c     &   mod_fourier(nloctyp)
 862 #else
 863       if (mod_tor) then
 864 #endif
 865       do i=-ntortyp+1,ntortyp-1
 866         do j=-ntortyp+1,ntortyp-1
 867           if (tor_mode.eq.0) then
 868           do iblock=1,2
 869             if (mask_tor(0,j,i,iblock).eq.1) then
 870               ii=ii+1
 871               weitor(0,j,i,iblock)=x(ii)
 872             else if (mask_tor(0,j,i,iblock).eq.2) then
 873               do k=1,nterm(j,i,iblock)
 874                 ii=ii+1
 875                 v1(k,j,i,iblock)=x(ii)
 876               enddo
 877             else if (mask_tor(0,j,i,iblock).gt.2) then
 878               do k=1,nterm(j,i,iblock)
 879                 ii=ii+1
 880                 v1(k,j,i,iblock)=x(ii)
 881               enddo
 882               do k=1,nterm(j,i,iblock)
 883                 ii=ii+1
 884                 v2(k,j,i,iblock)=x(ii)
 885               enddo
 886             endif
 887             v0ij=0.0d0
 888             si=-1.0d0
 889             do k=1,nterm(i,j,iblock)
 890               v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 891               v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 892               v0ij=v0ij+si*v1(k,i,j,iblock)
 893               si=-si
 894             enddo
 895             do k=1,nlor(i,j,iblock)
 896               v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 897             enddo
 898             v0(i,j,iblock)=v0ij
 899             v0(-i,-j,iblock)=v0ij
 900           enddo
 901           else if (tor_mode.eq.1) then
 902             if (mask_tor(0,j,i,1).eq.1) then
 903               ii=ii+1
 904               weitor(0,j,i,1)=x(ii)
 905             else if (mask_tor(0,j,i,1).eq.2) then
 906               do k=1,nterm_kcc(j,i)
 907                 do l=1,nterm_kcc_Tb(j,i)
 908                   do ll=1,nterm_kcc_Tb(j,i)
 909                     ii=ii+1
 910                     v1_kcc(ll,l,k,j,i)=x(ii)
 911                   enddo ! ll
 912                 enddo ! l
 913               enddo
 914             else if (mask_tor(0,j,i,1).gt.2) then
 915               do k=1,nterm_kcc(j,i)
 916                 do l=1,nterm_kcc_Tb(j,i)
 917                   do ll=1,nterm_kcc_Tb(j,i)
 918                     ii=ii+1
 919                     v1_kcc(ll,l,k,j,i)=x(ii)
 920                   enddo ! ll
 921                 enddo ! l
 922               enddo
 923               do k=1,nterm_kcc(j,i)
 924                 do l=1,nterm_kcc_Tb(j,i)
 925                   do ll=1,nterm_kcc_Tb(j,i)
 926                     ii=ii+1
 927                     v2_kcc(ll,l,k,j,i)=x(ii)
 928                   enddo! ll
 929                 enddo ! l
 930               enddo
 931             endif
 932           else
 933 #ifdef NEWCORR
 934 c Compute torsional coefficients from local energy surface expansion coeffs
 935             if (mask_tor(0,j,i,1).eq.1) then
 936               ii=ii+1
 937               weitor(0,j,i,1)=x(ii)
 938 c            else if (mask_tor(0,j,i,1).eq.2) then
 939             endif
 940 #else
 941             write (iout,*) "ERROR: TOR_MODE>1 ALLOWED ONLY WITH NEWCORR"
 942             stop
 943 #endif
 944           endif
 945         enddo
 946       enddo
 947
 948       endif
 949       if (mod_sccor) then
 950
 951       do i=-nsccortyp,nsccortyp
 952         do j=-nsccortyp,nsccortyp
 953           do k=1,3
 954           if (mask_tor(k,j,i,1).eq.1) then
 955             ii=ii+1
 956             weitor(k,j,i,1)=x(ii)
 957           else if (mask_tor(k,j,i,1).eq.2) then
 958             do l=1,nterm(j,i,1)
 959               ii=ii+1
 960               v1sccor(l,k,j,i)=x(ii)
 961             enddo
 962           else if (mask_tor(k,j,i,1).gt.2) then
 963             do l=1,nterm_sccor(j,i)
 964               ii=ii+1
 965               v1sccor(l,k,j,i)=x(ii)
 966             enddo
 967             do l=1,nterm_sccor(j,i)
 968               ii=ii+1
 969               v2sccor(l,k,j,i)=x(ii)
 970             enddo
 971           endif
 972           enddo
 973         enddo
 974       enddo
 975
 976       endif
 977
 978 c 7/8/17 AL Optimization of the bending potentials
 979
 980       if (mod_ang) then
 981         do i=0,nthetyp-1
 982           if (mask_ang(i).eq.0) cycle
 983           do j=1,nbend_kcc_TB(i)
 984             ii=ii+1
 985             v1bend_chyb(j,i)=x(ii)
 986             v1bend_chyb(j,-i)=x(ii)
 987           enddo
 988         enddo
 989       endif
 990
 991
 992       if (mod_fourier(nloctyp).gt.0) then
 993
 994 #ifdef NEWCORR
 995       do i=0,nloctyp-1
 996
 997         if (mod_fourier(i).gt.0) then
 998
 999           do j=1,2
1000             do k=1,3
1001               if (mask_bnew1(k,j,i).gt.0) then
1002                 ii=ii+1
1003                 bnew1(k,j,i)=x(ii)
1004               endif
1005             enddo
1006           enddo
1007
1008           do j=1,2
1009             do k=1,3
1010               if (mask_bnew2(k,j,i).gt.0) then
1011                 ii=ii+1
1012                 bnew2(k,j,i)=x(ii)
1013               endif
1014             enddo
1015           enddo
1016
1017           do j=1,2
1018             do k=1,3
1019               if (mask_ccnew(k,j,i).gt.0) then
1020                 ii=ii+1
1021                 ccnew(k,j,i)=x(ii)
1022               endif
1023             enddo
1024           enddo
1025
1026           do j=1,2
1027             do k=1,3
1028               if (mask_ddnew(k,j,i).gt.0) then
1029                 ii=ii+1
1030                 ddnew(k,j,i)=x(ii)
1031               endif
1032             enddo
1033           enddo
1034
1035           do j=1,3
1036             if (mask_e0new(j,i).gt.0) then
1037               ii=ii+1
1038               e0new(j,i)=x(ii)
1039             endif
1040           enddo
1041
1042           do j=1,2
1043             do k=1,2
1044               do l=1,2
1045                 if (mask_eenew(l,k,j,i).gt.0) then
1046                   ii=ii+1
1047                   eenew(l,k,j,i)=x(ii)
1048                 endif
1049               enddo
1050             enddo
1051           enddo
1052
1053         endif
1054
1055       enddo
1056
1057       do i=1,nloctyp-1
1058         do j=1,3
1059           bnew1(j,1,-i)= bnew1(j,1,i)
1060           bnew1(j,2,-i)=-bnew1(j,2,i)
1061           bnew2(j,1,-i)= bnew2(j,1,i)
1062           bnew2(j,2,-i)=-bnew2(j,2,i)
1063         enddo
1064         do j=1,3
1065 c          ccnew(j,1,i)=ccnew(j,1,i)/2
1066 c          ccnew(j,2,i)=ccnew(j,2,i)/2
1067 c          ddnew(j,1,i)=ddnew(j,1,i)/2
1068 c          ddnew(j,2,i)=ddnew(j,2,i)/2
1069           ccnew(j,1,-i)=ccnew(j,1,i)
1070           ccnew(j,2,-i)=-ccnew(j,2,i)
1071           ddnew(j,1,-i)=ddnew(j,1,i)
1072           ddnew(j,2,-i)=-ddnew(j,2,i)
1073         enddo
1074         e0new(1,-i)= e0new(1,i)
1075         e0new(2,-i)=-e0new(2,i)
1076         e0new(3,-i)=-e0new(3,i)
1077         do kk=1,2
1078           eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1079           eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1080           eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1081           eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1082         enddo
1083       enddo
1084       if (lprint) then
1085         write (iout,'(a)') "Coefficients of the multibody terms"
1086         do i=-nloctyp+1,nloctyp-1
1087           write (iout,*) "Type: ",onelet(iloctyp(i))
1088           write (iout,*) "Coefficients of the expansion of B1"
1089           do j=1,2
1090             write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1091           enddo
1092           write (iout,*) "Coefficients of the expansion of B2"
1093           do j=1,2
1094             write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1095           enddo
1096           write (iout,*) "Coefficients of the expansion of C"
1097           write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1098           write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1099           write (iout,*) "Coefficients of the expansion of D"
1100           write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1101           write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1102           write (iout,*) "Coefficients of the expansion of E"
1103           write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1104           do j=1,2
1105             do k=1,2
1106               write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1107             enddo
1108           enddo
1109         enddo
1110         call flush(iout)
1111       endif
1112
1113       IF (SPLIT_FOURIERTOR .and. tor_mode.eq.2) THEN
1114       do i=0,nloctyp-1
1115
1116         if (mod_fouriertor(i).gt.0) then
1117
1118           do j=1,2
1119             do k=1,3
1120               if (mask_bnew1tor(k,j,i).gt.0) then
1121                 ii=ii+1
1122                 bnew1tor(k,j,i)=x(ii)
1123               endif
1124             enddo
1125           enddo
1126
1127           do j=1,2
1128             do k=1,3
1129               if (mask_bnew2tor(k,j,i).gt.0) then
1130                 ii=ii+1
1131                 bnew2tor(k,j,i)=x(ii)
1132               endif
1133             enddo
1134           enddo
1135
1136           do j=1,2
1137             do k=1,3
1138               if (mask_ccnewtor(k,j,i).gt.0) then
1139                 ii=ii+1
1140                 ccnewtor(k,j,i)=x(ii)
1141               endif
1142             enddo
1143           enddo
1144
1145           do j=1,2
1146             do k=1,3
1147               if (mask_ddnewtor(k,j,i).gt.0) then
1148                 ii=ii+1
1149                 ddnewtor(k,j,i)=x(ii)
1150               endif
1151             enddo
1152           enddo
1153
1154           do j=1,3
1155             if (mask_e0newtor(j,i).gt.0) then
1156               ii=ii+1
1157               e0newtor(j,i)=x(ii)
1158             endif
1159           enddo
1160
1161           do j=1,2
1162             do k=1,2
1163               do l=1,2
1164                 if (mask_eenewtor(l,k,j,i).gt.0) then
1165                   ii=ii+1
1166                   eenewtor(l,k,j,i)=x(ii)
1167                 endif
1168               enddo
1169             enddo
1170           enddo
1171
1172         endif
1173
1174       enddo
1175
1176
1177       do i=1,nloctyp-1
1178         do j=1,3
1179           bnew1tor(j,1,-i)= bnew1tor(j,1,i)
1180           bnew1tor(j,2,-i)=-bnew1tor(j,2,i)
1181           bnew2tor(j,1,-i)= bnew2tor(j,1,i)
1182           bnew2tor(j,2,-i)=-bnew2tor(j,2,i)
1183         enddo
1184         do j=1,3
1185 c          ccnew(j,1,i)=ccnew(j,1,i)/2
1186 c          ccnew(j,2,i)=ccnew(j,2,i)/2
1187 c          ddnew(j,1,i)=ddnew(j,1,i)/2
1188 c          ddnew(j,2,i)=ddnew(j,2,i)/2
1189           ccnewtor(j,1,-i)=ccnewtor(j,1,i)
1190           ccnewtor(j,2,-i)=-ccnewtor(j,2,i)
1191           ddnewtor(j,1,-i)=ddnewtor(j,1,i)
1192           ddnewtor(j,2,-i)=-ddnewtor(j,2,i)
1193         enddo
1194         e0newtor(1,-i)= e0newtor(1,i)
1195         e0newtor(2,-i)=-e0newtor(2,i)
1196         e0newtor(3,-i)=-e0newtor(3,i)
1197         do kk=1,2
1198           eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1199           eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1200           eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1201           eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1202         enddo
1203       enddo
1204       if (lprint) then
1205         write (iout,'(a)')
1206      &     "Coefficients of the single-body torsional terms"
1207         do i=-nloctyp+1,nloctyp-1
1208           write (iout,*) "Type: ",onelet(iloctyp(i))
1209           write (iout,*) "Coefficients of the expansion of B1"
1210           do j=1,2
1211             write (iout,'(3hB1(,i1,1h),3f10.5)')
1212      &        j,(bnew1tor(k,j,i),k=1,3)
1213           enddo
1214           write (iout,*) "Coefficients of the expansion of B2"
1215           do j=1,2
1216             write (iout,'(3hB2(,i1,1h),3f10.5)')
1217      &        j,(bnew2tor(k,j,i),k=1,3)
1218           enddo
1219           write (iout,*) "Coefficients of the expansion of C"
1220           write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1221           write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1222           write (iout,*) "Coefficients of the expansion of D"
1223           write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1224           write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1225           write (iout,*) "Coefficients of the expansion of E"
1226           write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1227           do j=1,2
1228             do k=1,2
1229               write (iout,'(1hE,2i1,2f10.5)')
1230      &          j,k,(eenewtor(l,j,k,i),l=1,2)
1231             enddo
1232           enddo
1233         enddo
1234         call flush(iout)
1235       endif
1236
1237       ELSE
1238
1239       do i=-nloctyp+1,nloctyp-1
1240         do k=1,3
1241           do j=1,2
1242             bnew1tor(k,j,i)=bnew1(k,j,i)
1243           enddo
1244         enddo
1245         do k=1,3
1246           do j=1,2
1247             bnew2tor(k,j,i)=bnew2(k,j,i)
1248           enddo
1249         enddo
1250         do k=1,3
1251           ccnewtor(k,1,i)=ccnew(k,1,i)
1252           ccnewtor(k,2,i)=ccnew(k,2,i)
1253           ddnewtor(k,1,i)=ddnew(k,1,i)
1254           ddnewtor(k,2,i)=ddnew(k,2,i)
1255         enddo
1256
1257         do j=1,3
1258           e0newtor(j,i)=e0new(j,i)
1259         enddo
1260
1261         do j=1,2
1262           do k=1,2
1263             do l=1,2
1264               eenewtor(l,k,j,i)=eenew(l,k,j,i)
1265             enddo
1266           enddo
1267         enddo
1268
1269       enddo
1270
1271       ENDIF ! SPLIT_FOURIERTOR
1272
1273       if (tor_mode.eq.2) then
1274 c Recalculate new torsional potential coefficients
1275         do i=-ntortyp+1,ntortyp-1
1276           do j=-ntortyp+1,ntortyp-1
1277             do k=1,nterm_kcc_Tb(j,i)
1278               do l=1,nterm_kcc_Tb(j,i)
1279                 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1280      &                         +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1281                 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1282      &                         +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1283               enddo
1284             enddo
1285             do k=1,nterm_kcc_Tb(j,i)
1286               do l=1,nterm_kcc_Tb(j,i)
1287 c                v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1288 c     &                         -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1289 c                v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1290 c     &                         +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1291                 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1292      &                         -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1293                 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1294      &                         +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1295               enddo
1296             enddo
1297           enddo
1298         enddo
1299       endif
1300 #else
1301       do i=0,nloctyp-1
1302
1303         if (mod_fourier(i).gt.0) then
1304
1305           do j=2,13
1306             if (mask_fourier(j,i).gt.0) then
1307               ii=ii+1
1308               b(j,i)=x(ii)
1309             endif
1310           enddo
1311
1312           if (lprint) then
1313             write (iout,*) "Fourier coefficients of cumulants"
1314             write (iout,*) 'Type',i
1315             write (iout,'(a,i2,a,f10.5)') ('b(',j,')=',b(j,i),j=1,13)
1316           endif
1317
1318           if (i.gt.0) then
1319             b(2,-i)= b(2,i)
1320             b(3,-i)= b(3,i)
1321             b(4,-i)=-b(4,i)
1322             b(5,-i)=-b(5,i)
1323           endif
1324 c          B1(1,i)  = b(3,i)
1325 c          B1(2,i)  = b(5,i)
1326 c          B1tilde(1,i) = b(3,i)
1327 c          B1tilde(2,i) =-b(5,i)
1328 c          B2(1,i)  = b(2,i)
1329 c          B2(2,i)  = b(4,i)
1330           CCold(1,1,i)= b(7,i)
1331           CCold(2,2,i)=-b(7,i)
1332           CCold(2,1,i)= b(9,i)
1333           CCold(1,2,i)= b(9,i)
1334           CCold(1,1,-i)= b(7,i)
1335           CCold(2,2,-i)=-b(7,i)
1336           CCold(2,1,-i)=-b(9,i)
1337           CCold(1,2,-i)=-b(9,i)
1338 c          Ctilde(1,1,i)=b(7,i)
1339 c          Ctilde(1,2,i)=b(9,i)
1340 c          Ctilde(2,1,i)=-b(9,i)
1341 c          Ctilde(2,2,i)=b(7,i)
1342           DDold(1,1,i)= b(6,i)
1343           DDold(2,2,i)=-b(6,i)
1344           DDold(2,1,i)= b(8,i)
1345           DDold(1,2,i)= b(8,i)
1346           DDold(1,1,-i)= b(6,i)
1347           DDold(2,2,-i)=-b(6,i)
1348           DDold(2,1,-i)=-b(8,i)
1349           DDold(1,2,-i)=-b(8,i)
1350 c          Dtilde(1,1,i)=b(6,i)
1351 c          Dtilde(1,2,i)=b(8,i)
1352 c          Dtilde(2,1,i)=-b(8,i)
1353 c          Dtilde(2,2,i)=b(6,i)
1354           EEold(1,1,i)= b(10,i)+b(11,i)
1355           EEold(2,2,i)=-b(10,i)+b(11,i)
1356           EEold(2,1,i)= b(12,i)-b(13,i)
1357           EEold(1,2,i)= b(12,i)+b(13,i)
1358           EEold(1,1,i)= b(10,i)+b(11,i)
1359           EEold(2,2,i)=-b(10,i)+b(11,i)
1360           EEold(2,1,i)= b(12,i)-b(13,i)
1361           EEold(1,2,i)= b(12,i)+b(13,i)
1362           EEold(1,1,-i)= b(10,i)+b(11,i)
1363           EEold(2,2,-i)=-b(10,i)+b(11,i)
1364           EEold(2,1,-i)=-b(12,i)+b(13,i)
1365           EEold(1,2,-i)=-b(12,i)-b(13,i)
1366
1367         endif
1368
1369       enddo
1370       if (lprint) then
1371       write (iout,*)
1372       write (iout,*)
1373      &"Coefficients of the cumulants (independent of valence angles)"
1374       do i=-nloctyp+1,nloctyp-1
1375         write (iout,*) 'Type ',onelet(iloctyp(i))
1376         write (iout,*) 'B1'
1377         write(iout,'(2f10.5)') B(3,i),B(5,i)
1378         write (iout,*) 'B2'
1379         write(iout,'(2f10.5)') B(2,i),B(4,i)
1380         write (iout,*) 'CC'
1381         do j=1,2
1382           write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1383         enddo
1384         write(iout,*) 'DD'
1385         do j=1,2
1386           write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1387         enddo
1388         write(iout,*) 'EE'
1389         do j=1,2
1390           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1391         enddo
1392       enddo
1393       endif
1394 #endif
1395
1396       endif
1397 c SC sigmas and anisotropies as parameters
1398       if (mod_side_other) then
1399 c        do i=1,ntyp
1400 c          write (iout,*) "mask",i,(mask_sigma(k,i),k=1,4)
1401 c        enddo
1402         do i=1,ntyp
1403           if (mask_sigma(1,i).gt.0) then
1404             ii=ii+1
1405             sigma0(i)=x(ii)
1406 c            write(iout,*) "sig0",i,ii,x(ii)
1407           endif
1408           if (mask_sigma(2,i).gt.0) then
1409             ii=ii+1
1410             sigii(i)=x(ii)
1411 c            write(iout,*) "sigi",i,ii,x(ii)
1412           endif
1413           if (mask_sigma(3,i).gt.0) then
1414             ii=ii+1
1415             chip0(i)=x(ii)
1416 c            write(iout,*) "chip",i,ii,x(ii)
1417           endif
1418           if (mask_sigma(4,i).gt.0) then
1419             ii=ii+1
1420             alp(i)=x(ii)
1421 c            write(iout,*) "alp",i,ii,x(ii)
1422           endif
1423         enddo
1424       endif
1425
1426       if (mod_side.or.mod_side_other) then
1427         if (lprint) write(iout,'(/3a,2i3)') 'Potential is ',
1428      &    potname(ipot),', exponents are ',expon,2*expon
1429 C-----------------------------------------------------------------------
1430 C Calculate the "working" parameters of SC interactions.
1431         dwa16=2.0d0**(1.0d0/6.0d0)
1432         if (ipot.eq.4) then
1433           do i=1,ntyp
1434             chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1435           enddo
1436         endif
1437         do i=1,ntyp
1438           do j=i,ntyp
1439             sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1440             sigma(j,i)=sigma(i,j)
1441             rs0(i,j)=dwa16*sigma(i,j)
1442             rs0(j,i)=rs0(i,j)
1443           enddo
1444         enddo
1445         if (lprint) call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1446         if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1447      &   'Working parameters of the SC interactions:',
1448      &   '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1449      &   '  chi1   ','   chi2   '
1450         do i=1,ntyp
1451           do j=i,ntyp
1452             epsij=eps(i,j)
1453             if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1454               rrij=sigma(i,j)
1455             else
1456               rrij=rr0(i)+rr0(j)
1457             endif
1458             r0(i,j)=rrij
1459             r0(j,i)=rrij
1460             rrij=rrij**expon
1461             epsij=eps(i,j)
1462 c           sigeps=dsign(1.0D0,epsij)
1463             aa(i,j)=ftune_eps(epsij)*rrij*rrij
1464             bb(i,j)=-epsij*rrij
1465             aa(j,i)=aa(i,j)
1466             bb(j,i)=bb(i,j)
1467             if (ipot.gt.2) then
1468               sigt1sq=sigma0(i)**2
1469               sigt2sq=sigma0(j)**2
1470               sigii1=sigii(i)
1471               sigii2=sigii(j)
1472               ratsig1=sigt2sq/sigt1sq
1473               ratsig2=1.0D0/ratsig1
1474               chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1475               if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1476               rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1477             else
1478               rsum_max=sigma(i,j)
1479             endif
1480             sigmaii(i,j)=rsum_max
1481             sigmaii(j,i)=rsum_max
1482             if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max)
1483      &      then
1484               r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1485               augm(i,j)=epsij*r_augm**(2*expon)
1486 c             augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1487               augm(j,i)=augm(i,j)
1488             else
1489               augm(i,j)=0.0D0
1490               augm(j,i)=0.0D0
1491             endif
1492             if (lprint) then
1493               write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1494      &        restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1495      &        sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1496             endif
1497           enddo
1498         enddo
1499       endif
1500
1501
1502       if (mod_elec) then
1503
1504       do i=1,2
1505         do j=1,i
1506           if (mask_elec(i,j,1).gt.0) then
1507             ii=ii+1
1508             epp(i,j)=x(ii)
1509             if (epp(i,j).lt.0.0d0) epp(i,j)=0.0d0
1510             epp(j,i)=x(ii)
1511           endif
1512           if (mask_elec(i,j,2).gt.0) then
1513             ii=ii+1
1514             rpp(i,j)=dabs(x(ii))
1515             rpp(j,i)=dabs(x(ii))
1516           endif
1517           if (mask_elec(i,j,3).gt.0) then
1518             ii=ii+1
1519             elpp6(i,j)=x(ii)
1520             if (elpp6(i,j).gt.0.0d0) elpp6(i,j)=0.0d0
1521             elpp6(j,i)=elpp6(i,j)
1522           endif
1523           if (mask_elec(i,j,4).gt.0) then
1524             ii=ii+1
1525             elpp3(i,j)=x(ii)
1526             elpp3(j,i)=x(ii)
1527           endif
1528           rri=rpp(i,j)**6
1529           app (i,j)=epp(i,j)*rri*rri
1530           app (j,i)=epp(j,i)*rri*rri
1531           bpp (i,j)=-2.0D0*epp(i,j)*rri
1532           bpp (j,i)=-2.0D0*epp(i,j)*rri
1533           ael6(i,j)=elpp6(i,j)*4.2D0**6
1534           ael6(j,i)=elpp6(i,j)*4.2D0**6
1535           ael3(i,j)=elpp3(i,j)*4.2D0**3
1536           ael3(j,i)=elpp3(i,j)*4.2D0**3
1537         enddo
1538       enddo
1539       if (lprint) then
1540         write (iout,'(/a)') 'Electrostatic interaction constants:'
1541         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1542      &            'IT','JT','APP','BPP','AEL6','AEL3'
1543         do i=1,2
1544           do j=1,2
1545             write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1546      &                    ael6(i,j),ael3(i,j)
1547           enddo
1548         enddo
1549       endif
1550
1551       endif
1552 C
1553 C Define the SC-p interaction constants
1554 C
1555       if (mod_scp) then
1556
1557       if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+
1558      &        mask_scp(0,2,2) .gt. 0) then
1559         do j=1,2
1560           if (mask_scp(0,j,1).gt.0) then
1561             ii=ii+1
1562             do i=1,20
1563               eps_scp(i,j)=x(ii)
1564             enddo
1565           endif
1566           if (mask_scp(0,j,2).gt.0) then
1567             ii=ii+1
1568             do i=1,20
1569               rscp(i,j)=x(ii)
1570             enddo
1571           endif
1572         enddo
1573       else
1574         do i=1,20
1575           do j=1,2
1576             if (mask_scp(i,j,1).gt.0) then
1577               ii=ii+1
1578               eps_scp(i,j)=x(ii)
1579             endif
1580             if (mask_scp(i,j,2).gt.0) then
1581               ii=ii+1
1582               rscp(i,j)=x(ii)
1583             endif
1584           enddo
1585         enddo
1586       endif
1587       do i=1,20
1588         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1589         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1590         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1591         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1592       enddo
1593
1594       if (lprint) then
1595         write (iout,*) "Parameters of SC-p interactions:"
1596         do i=1,20
1597           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1598      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1599         enddo
1600       endif
1601
1602       endif
1603
1604 c Weight of the shielding term
1605       if (imask(ind_shield).eq.1) then
1606         ii=ii+1
1607         ww(ind_shield)=x(ii)
1608       endif
1609       wshield=ww(ind_shield)
1610
1611 c      write (iout,*) "x2w: nvar",nvarr
1612
1613 c Base of heat capacity
1614
1615       do iprot=1,nprot
1616         if (nbeta(2,iprot).gt.0) then
1617           ii=ii+1
1618           heatbase(iprot)=x(ii)
1619         endif
1620       enddo
1621
1622       if (ii.ne.nvarr) then
1623         write (iout,*) "!!!!! SEVERE ERROR IN X2W !!!!! (chuj nastapil)"
1624         write (iout,*) "NVAR=",nvarr," LAST VARIABLE INDEX=",ii
1625         stop "CHUJ ABSOLUTNY!!!!"
1626       endif
1627
1628       return
1629       end
1630 c---------------------------------------------------------------------------
1631       double precision function ftune_eps(x)
1632       implicit none
1633       double precision x,y
1634       double precision delta /1.0d-5/
1635       if (dabs(x).lt.delta) then
1636         y=x/delta
1637         ftune_eps=delta*(3.0d0+6.0d0*y**2-y**4)/8
1638       else
1639         ftune_eps=dabs(x)
1640       endif
1641       return
1642       end
1643 c---------------------------------------------------------------------------
1644       double precision function ftune_epsprim(x)
1645       implicit none
1646       double precision x,y
1647       double precision delta /1.0d-5/
1648       if (dabs(x).lt.delta) then
1649         y=x/delta
1650         ftune_epsprim=1.5d0*y-0.5*y**3
1651       else
1652         ftune_epsprim=dsign(1.0d0,x)
1653       endif
1654       return
1655       end
1656 c---------------------------------------------------------------------------
1657       double precision function ftune_epsbis(x)
1658       implicit none
1659       double precision x,y
1660       double precision delta /1.0d-5/
1661       if (dabs(x).lt.delta) then
1662         y=x/delta
1663         ftune_epsbis=(1.5d0-1.5d0*y**3)/delta
1664       else
1665         ftune_epsbis=0.0d0
1666       endif
1667       return
1668       end
1669 c---------------------------------------------------------------------------
1670       logical function in_sc_pair_list(iti,itj)
1671       implicit none
1672       include "DIMENSIONS"
1673       include "DIMENSIONS.ZSCOPT"
1674       include "COMMON.WEIGHTS"
1675       integer i,iti,itj
1676       do i=1,npair_sc
1677         if (ityp_psc(1,i).eq.iti .and. ityp_psc(2,i).eq.itj .or.
1678      &      ityp_psc(1,i).eq.itj .and. ityp_psc(2,i).eq.iti) then
1679           in_sc_pair_list=.true.
1680           return
1681         endif
1682       enddo
1683       in_sc_pair_list=.false.
1684       return
1685       end