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