source/unres/src_MD-M-SAXS-homology/readrtns_CSA.F

   1       subroutine readrtns
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.SETUP'
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.SBRIDGE'
  10       include 'COMMON.IOUNITS'
  11       include 'COMMON.SPLITELE'
  12       logical file_exist
  13 C Read job setup parameters
  14       call read_control
  15 C Read force-field parameters except weights
  16       call parmread
  17 C Read control parameters for energy minimzation if required
  18       if (minim) call read_minim
  19 C Read MCM control parameters if required
  20       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  21 C Read MD control parameters if reqjuired
  22       if (modecalc.eq.12) call read_MDpar
  23 C Read MREMD control parameters if required
  24       if (modecalc.eq.14) then
  25          call read_MDpar
  26          call read_REMDpar
  27       endif
  28 C Read MUCA control parameters if required
  29       if (lmuca) call read_muca
  30 C Read CSA control parameters if required (from fort.40 if exists
  31 C otherwise from general input file)
  32       if (modecalc.eq.8) then
  33        inquire (file="fort.40",exist=file_exist)
  34        if (.not.file_exist) call csaread
  35       endif
  36 cfmc      if (modecalc.eq.10) call mcmfread
  37 c Read energy-term weights and disulfide parameters
  38       call weightread
  39 C Read molecule information, molecule geometry, energy-term weights, and
  40 C restraints if requested
  41       call molread
  42 C Print restraint information
  43 #ifdef MPI
  44       if (.not. out1file .or. me.eq.king) then
  45 #endif
  46       if (nhpb.gt.nss)  then
  47         write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
  48      &  "The following",nhpb-nss,
  49      &  " distance restraints have been imposed:",
  50      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
  51      &  "  score"," type"
  52         do i=nss+1,nhpb
  53           write (iout,'(4i5,2f8.2,3f10.5,2i5)')i-nss,ihpb(i),jhpb(i),
  54      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
  55      &  irestr_type(i)
  56         enddo
  57       endif
  58       if (npeak.gt.0) then
  59         write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
  60      &  "The following",npeak,
  61      &  " NMR peak restraints have been imposed:",
  62      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
  63      &  "  score"," type"," ipeak"
  64         do i=1,npeak
  65           do j=ipeak(1,i),ipeak(2,i)
  66             write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
  67      &       jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
  68      &       forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
  69           enddo
  70         enddo
  71         write (iout,*) "The ipeak array"
  72         do i=1,npeak
  73           write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
  74         enddo
  75       endif
  76 #ifdef MPI
  77       endif
  78 #endif
  79 c      print *,"Processor",myrank," leaves READRTNS"
  80       return
  81       end
  82 C-------------------------------------------------------------------------------
  83       subroutine read_control
  84 C
  85 C Read contorl data
  86 C
  87       implicit real*8 (a-h,o-z)
  88       include 'DIMENSIONS'
  89 #ifdef MP
  90       include 'mpif.h'
  91       logical OKRandom, prng_restart
  92       real*8  r1
  93 #endif
  94       include 'COMMON.IOUNITS'
  95       include 'COMMON.TIME1'
  96       include 'COMMON.THREAD'
  97       include 'COMMON.SBRIDGE'
  98       include 'COMMON.CONTROL'
  99       include 'COMMON.MCM'
 100       include 'COMMON.MAP'
 101       include 'COMMON.HEADER'
 102       include 'COMMON.CSA'
 103       include 'COMMON.CHAIN'
 104       include 'COMMON.MUCA'
 105       include 'COMMON.MD'
 106       include 'COMMON.FFIELD'
 107       include 'COMMON.INTERACT'
 108       include 'COMMON.SETUP'
 109       include 'COMMON.SPLITELE'
 110       include 'COMMON.SHIELD'
 111       include 'COMMON.GEO'
 112       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
 113       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
 114       character*80 ucase
 115       character*320 controlcard
 116
 117       nglob_csa=0
 118       eglob_csa=1d99
 119       nmin_csa=0
 120       read (INP,'(a)') titel
 121       call card_concat(controlcard)
 122 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
 123 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
 124       call reada(controlcard,'SEED',seed,0.0D0)
 125       call random_init(seed)
 126 C Set up the time limit (caution! The time must be input in minutes!)
 127       read_cart=index(controlcard,'READ_CART').gt.0
 128       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
 129 C this variable with_theta_constr is the variable which allow to read and execute the
 130 C constrains on theta angles WITH_THETA_CONSTR is the keyword
 131       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
 132       write (iout,*) "with_theta_constr ",with_theta_constr
 133       call readi(controlcard,'NSAXS',nsaxs,0)
 134       call readi(controlcard,'SAXS_MODE',saxs_mode,0)
 135       call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
 136       call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
 137       write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
 138      &   SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
 139       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
 140       read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
 141       out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
 142       out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
 143       call readi(controlcard,'SYM',symetr,1)
 144       call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
 145       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 146       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 147       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 148       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 149       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 150       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 151       call reada(controlcard,'DRMS',drms,0.1D0)
 152       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 153        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 154        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 155        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 156        write (iout,'(a,f10.1)')'DRMS    = ',drms
 157        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 158        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 159       endif
 160       call readi(controlcard,'NZ_START',nz_start,0)
 161       call readi(controlcard,'NZ_END',nz_end,0)
 162 c      call readi(controlcard,'IZ_SC',iz_sc,0)
 163       timlim=60.0D0*timlim
 164       safety = 60.0d0*safety
 165       timem=timlim
 166       modecalc=0
 167       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 168       minim=(index(controlcard,'MINIMIZE').gt.0)
 169       dccart=(index(controlcard,'CART').gt.0)
 170       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 171       overlapsc=.not.overlapsc
 172       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 173       searchsc=.not.searchsc
 174       sideadd=(index(controlcard,'SIDEADD').gt.0)
 175       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 176       mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
 177       outpdb=(index(controlcard,'PDBOUT').gt.0)
 178       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 179       pdbref=(index(controlcard,'PDBREF').gt.0)
 180       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 181       indpdb=index(controlcard,'PDBSTART')
 182       extconf=(index(controlcard,'EXTCONF').gt.0)
 183       AFMlog=(index(controlcard,'AFM'))
 184       selfguide=(index(controlcard,'SELFGUIDE'))
 185 c      print *,'AFMlog',AFMlog,selfguide,"KUPA"
 186       call readi(controlcard,'TUBEMOD',tubelog,0)
 187 c      write (iout,*) TUBElog,"TUBEMODE"
 188       call readi(controlcard,'IPRINT',iprint,0)
 189 C SHIELD keyword sets if the shielding effect of side-chains is used
 190 C 0 denots no shielding is used all peptide are equally despite the
 191 C solvent accesible area
 192 C 1 the newly introduced function
 193 C 2 reseved for further possible developement
 194       call readi(controlcard,'SHIELD',shield_mode,0)
 195 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 196         write(iout,*) "shield_mode",shield_mode
 197 C      endif
 198       call readi(controlcard,'TORMODE',tor_mode,0)
 199 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 200         write(iout,*) "torsional and valence angle mode",tor_mode
 201       call readi(controlcard,'MAXGEN',maxgen,10000)
 202       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 203       call readi(controlcard,"KDIAG",kdiag,0)
 204       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 205       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 206      & write (iout,*) "RESCALE_MODE",rescale_mode
 207       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 208       if (index(controlcard,'REGULAR').gt.0.0D0) then
 209         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 210         modecalc=1
 211         refstr=.true.
 212       endif
 213       if (refstr) then
 214       if (index(controlcard,"CASC").gt.0) then
 215         iz_sc=1
 216 c      else if (index(controlcard,"CAONLY").gt.0) then
 217 c        iz_sc=0
 218       else if (index(controlcard,"SCONLY").gt.0) then
 219         iz_sc=2
 220       else
 221         iz_sc=0
 222 c        write (iout,*) "ERROR! Must specify CASC CAONLY or SCONLY when",
 223 c     &    " specifying REFSTR, PDBREF or REGULAR."
 224 c        stop
 225       endif
 226       endif
 227       if (index(controlcard,'CHECKGRAD').gt.0) then
 228         modecalc=5
 229         if (index(controlcard,' CART').gt.0) then
 230           icheckgrad=1
 231         elseif (index(controlcard,'CARINT').gt.0) then
 232           icheckgrad=2
 233         else
 234           icheckgrad=3
 235         endif
 236         call reada(controlcard,'DELTA',aincr,1.0d-7)
 237 c        write (iout,*) "icheckgrad",icheckgrad
 238       elseif (index(controlcard,'THREAD').gt.0) then
 239         modecalc=2
 240         call readi(controlcard,'THREAD',nthread,0)
 241         if (nthread.gt.0) then
 242           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 243         else
 244           if (fg_rank.eq.0)
 245      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 246           stop 'Error termination in Read_Control.'
 247         endif
 248       else if (index(controlcard,'MCMA').gt.0) then
 249         modecalc=3
 250       else if (index(controlcard,'MCEE').gt.0) then
 251         modecalc=6
 252       else if (index(controlcard,'MULTCONF').gt.0) then
 253         modecalc=4
 254       else if (index(controlcard,'MAP').gt.0) then
 255         modecalc=7
 256         call readi(controlcard,'MAP',nmap,0)
 257       else if (index(controlcard,'CSA').gt.0) then
 258         modecalc=8
 259 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 260 crc
 261 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 262 crc
 263 crc        modecalc=9
 264 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 265 cfmc        modecalc=10
 266       else if (index(controlcard,'SOFTREG').gt.0) then
 267         modecalc=11
 268       else if (index(controlcard,'CHECK_BOND').gt.0) then
 269         modecalc=-1
 270       else if (index(controlcard,'TEST').gt.0) then
 271         modecalc=-2
 272       else if (index(controlcard,'MD').gt.0) then
 273         modecalc=12
 274       else if (index(controlcard,'RE ').gt.0) then
 275         modecalc=14
 276       endif
 277
 278       lmuca=index(controlcard,'MUCA').gt.0
 279       call readi(controlcard,'MUCADYN',mucadyn,0)
 280       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 281       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 282      & then
 283        write (iout,*) 'MUCADYN=',mucadyn
 284        write (iout,*) 'MUCASMOOTH=',muca_smooth
 285       endif
 286
 287       iscode=index(controlcard,'ONE_LETTER')
 288       indphi=index(controlcard,'PHI')
 289       indback=index(controlcard,'BACK')
 290       iranconf=index(controlcard,'RAND_CONF')
 291       i2ndstr=index(controlcard,'USE_SEC_PRED')
 292       gradout=index(controlcard,'GRADOUT').gt.0
 293       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 294 C  DISTCHAINMAX become obsolete for periodic boundry condition
 295       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
 296 C Reading the dimensions of box in x,y,z coordinates
 297       call reada(controlcard,'BOXX',boxxsize,100.0d0)
 298       call reada(controlcard,'BOXY',boxysize,100.0d0)
 299       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
 300 c Cutoff range for interactions
 301       call reada(controlcard,"R_CUT",r_cut,15.0d0)
 302       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 303       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
 304       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
 305       if (lipthick.gt.0.0d0) then
 306        bordliptop=(boxzsize+lipthick)/2.0
 307        bordlipbot=bordliptop-lipthick
 308 C      endif
 309       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
 310      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
 311       buflipbot=bordlipbot+lipbufthick
 312       bufliptop=bordliptop-lipbufthick
 313       if ((lipbufthick*2.0d0).gt.lipthick)
 314      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
 315       endif
 316       write(iout,*) "bordliptop=",bordliptop
 317       write(iout,*) "bordlipbot=",bordlipbot
 318       write(iout,*) "bufliptop=",bufliptop
 319       write(iout,*) "buflipbot=",buflipbot
 320       write (iout,*) "SHIELD MODE",shield_mode
 321       if (TUBElog.gt.0) then
 322        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
 323        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
 324        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
 325        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
 326        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
 327        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
 328        buftubebot=bordtubebot+tubebufthick
 329        buftubetop=bordtubetop-tubebufthick
 330       endif
 331 c      if (shield_mode.gt.0) then
 332 c      pi=3.141592d0
 333 C VSolvSphere the volume of solving sphere
 334 C      print *,pi,"pi"
 335 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
 336 C there will be no distinction between proline peptide group and normal peptide
 337 C group in case of shielding parameters
 338 c      write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
 339 c      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
 340 c      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
 341 c      write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
 342 c     &  VSolvSphere_div
 343 C long axis of side chain
 344 c      do i=1,ntyp
 345 c      long_r_sidechain(i)=vbldsc0(1,i)
 346 c      short_r_sidechain(i)=sigma0(i)
 347 c      enddo
 348 c      buff_shield=1.0d0
 349 c      endif
 350       if (me.eq.king .or. .not.out1file )
 351      &  write (iout,*) "DISTCHAINMAX",distchainmax
 352
 353       if(me.eq.king.or..not.out1file)
 354      & write (iout,'(2a)') diagmeth(kdiag),
 355      &  ' routine used to diagonalize matrices.'
 356       return
 357       end
 358 c--------------------------------------------------------------------------
 359       subroutine read_REMDpar
 360 C
 361 C Read REMD settings
 362 C
 363       implicit real*8 (a-h,o-z)
 364       include 'DIMENSIONS'
 365       include 'COMMON.IOUNITS'
 366       include 'COMMON.TIME1'
 367       include 'COMMON.MD'
 368 #ifndef LANG0
 369       include 'COMMON.LANGEVIN'
 370 #else
 371       include 'COMMON.LANGEVIN.lang0'
 372 #endif
 373       include 'COMMON.INTERACT'
 374       include 'COMMON.NAMES'
 375       include 'COMMON.GEO'
 376       include 'COMMON.REMD'
 377       include 'COMMON.CONTROL'
 378       include 'COMMON.SETUP'
 379       character*80 ucase
 380       character*320 controlcard
 381       character*3200 controlcard1
 382       integer iremd_m_total
 383
 384       if(me.eq.king.or..not.out1file)
 385      & write (iout,*) "REMD setup"
 386
 387       call card_concat(controlcard)
 388       call readi(controlcard,"NREP",nrep,3)
 389       call readi(controlcard,"NSTEX",nstex,1000)
 390       call reada(controlcard,"RETMIN",retmin,10.0d0)
 391       call reada(controlcard,"RETMAX",retmax,1000.0d0)
 392       mremdsync=(index(controlcard,'SYNC').gt.0)
 393       call readi(controlcard,"NSYN",i_sync_step,100)
 394       restart1file=(index(controlcard,'REST1FILE').gt.0)
 395       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
 396       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
 397       if(max_cache_traj_use.gt.max_cache_traj)
 398      &           max_cache_traj_use=max_cache_traj
 399       if(me.eq.king.or..not.out1file) then
 400 cd       if (traj1file) then
 401 crc caching is in testing - NTWX is not ignored
 402 cd        write (iout,*) "NTWX value is ignored"
 403 cd        write (iout,*) "  trajectory is stored to one file by master"
 404 cd        write (iout,*) "  before exchange at NSTEX intervals"
 405 cd       endif
 406        write (iout,*) "NREP= ",nrep
 407        write (iout,*) "NSTEX= ",nstex
 408        write (iout,*) "SYNC= ",mremdsync
 409        write (iout,*) "NSYN= ",i_sync_step
 410        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
 411       endif
 412       remd_tlist=.false.
 413       if (index(controlcard,'TLIST').gt.0) then
 414          remd_tlist=.true.
 415          call card_concat(controlcard1)
 416          read(controlcard1,*) (remd_t(i),i=1,nrep)
 417          if(me.eq.king.or..not.out1file)
 418      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep)
 419       endif
 420       remd_mlist=.false.
 421       if (index(controlcard,'MLIST').gt.0) then
 422          remd_mlist=.true.
 423          call card_concat(controlcard1)
 424          read(controlcard1,*) (remd_m(i),i=1,nrep)
 425          if(me.eq.king.or..not.out1file) then
 426           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
 427           iremd_m_total=0
 428           do i=1,nrep
 429            iremd_m_total=iremd_m_total+remd_m(i)
 430           enddo
 431            write (iout,*) 'Total number of replicas ',iremd_m_total
 432           endif
 433          endif
 434         if (adaptive) then
 435           write (iout,*)
 436      &      "Adaptive (PMF-biased) umbrella sampling will be run"
 437           call PMFread
 438         endif
 439       if(me.eq.king.or..not.out1file)
 440      &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
 441       return
 442       end
 443 c--------------------------------------------------------------------------
 444       subroutine read_MDpar
 445 C
 446 C Read MD settings
 447 C
 448       implicit real*8 (a-h,o-z)
 449       include 'DIMENSIONS'
 450       include 'COMMON.IOUNITS'
 451       include 'COMMON.TIME1'
 452       include 'COMMON.MD'
 453 #ifndef LANG0
 454       include 'COMMON.LANGEVIN'
 455 #else
 456       include 'COMMON.LANGEVIN.lang0'
 457 #endif
 458       include 'COMMON.INTERACT'
 459       include 'COMMON.NAMES'
 460       include 'COMMON.GEO'
 461       include 'COMMON.SETUP'
 462       include 'COMMON.CONTROL'
 463       include 'COMMON.SPLITELE'
 464       include 'COMMON.FFIELD'
 465       character*80 ucase
 466       character*320 controlcard
 467
 468       call card_concat(controlcard)
 469       call readi(controlcard,"NSTEP",n_timestep,1000000)
 470       call readi(controlcard,"NTWE",ntwe,100)
 471       call readi(controlcard,"NTWX",ntwx,1000)
 472       call reada(controlcard,"DT",d_time,1.0d-1)
 473       call reada(controlcard,"DVMAX",dvmax,2.0d1)
 474       call reada(controlcard,"DAMAX",damax,1.0d1)
 475       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
 476       call readi(controlcard,"LANG",lang,0)
 477       RESPA = index(controlcard,"RESPA") .gt. 0
 478       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
 479       ntime_split0=ntime_split
 480       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
 481       ntime_split0=ntime_split
 482 c      call reada(controlcard,"R_CUT",r_cut,2.0d0)
 483 c      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 484       rest = index(controlcard,"REST").gt.0
 485       tbf = index(controlcard,"TBF").gt.0
 486       usampl = index(controlcard,"USAMPL").gt.0
 487       scale_umb =  index(controlcard,"SCALE_UMB").gt.0
 488       adaptive = index(controlcard,"ADAPTIVE").gt.0
 489       mdpdb = index(controlcard,"MDPDB").gt.0
 490       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 491       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
 492       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
 493       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
 494       if (count_reset_moment.eq.0) count_reset_moment=1000000000
 495       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
 496       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
 497       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
 498       if (count_reset_vel.eq.0) count_reset_vel=1000000000
 499       large = index(controlcard,"LARGE").gt.0
 500       print_compon = index(controlcard,"PRINT_COMPON").gt.0
 501       rattle = index(controlcard,"RATTLE").gt.0
 502       preminim = index(controlcard,"PREMINIM").gt.0
 503       if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then
 504         preminim=.true.
 505         dccart=.true.
 506       endif
 507       write (iout,*) "PREMINIM ",preminim
 508       if (preminim) then
 509         if (index(controlcard,'CART').gt.0) dccart=.true.
 510         write (iout,*) "dccart ",dccart
 511         write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR')
 512         if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim
 513       endif
 514 c  if performing umbrella sampling, fragments constrained are read from the fragment file
 515       nset=0
 516       if(usampl) then
 517         write (iout,*) "Umbrella sampling will be run"
 518         if (scale_umb.and.adaptive) then
 519           write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
 520           write (iout,*) "Select one of those and re-run the job."
 521           stop
 522         endif
 523         if (scale_umb) write (iout,*)
 524      &"Umbrella-restraint force constants will be scaled by temperature"
 525         call read_fragments
 526       endif
 527
 528       if(me.eq.king.or..not.out1file) then
 529        write (iout,*)
 530        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
 531        write (iout,*)
 532        write (iout,'(a)') "The units are:"
 533        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
 534        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
 535      &  " acceleration: angstrom/(48.9 fs)**2"
 536        write (iout,'(a)') "energy: kcal/mol, temperature: K"
 537        write (iout,*)
 538        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
 539        write (iout,'(a60,f10.5,a)')
 540      &  "Initial time step of numerical integration:",d_time,
 541      &  " natural units"
 542        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
 543        if (RESPA) then
 544         write (iout,'(2a,i4,a)')
 545      &    "A-MTS algorithm used; initial time step for fast-varying",
 546      &    " short-range forces split into",ntime_split," steps."
 547         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
 548      &   r_cut," lambda",rlamb
 549        endif
 550        write (iout,'(2a,f10.5)')
 551      &  "Maximum acceleration threshold to reduce the time step",
 552      &  "/increase split number:",damax
 553        write (iout,'(2a,f10.5)')
 554      &  "Maximum predicted energy drift to reduce the timestep",
 555      &  "/increase split number:",edriftmax
 556        write (iout,'(a60,f10.5)')
 557      & "Maximum velocity threshold to reduce velocities:",dvmax
 558        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
 559        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
 560        if (rattle) write (iout,'(a60)')
 561      &  "Rattle algorithm used to constrain the virtual bonds"
 562       endif
 563       reset_fricmat=1000
 564       if (lang.gt.0) then
 565         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
 566         call reada(controlcard,"RWAT",rwat,1.4d0)
 567         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
 568         surfarea=index(controlcard,"SURFAREA").gt.0
 569         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
 570         if(me.eq.king.or..not.out1file)then
 571          write (iout,'(/a,$)') "Langevin dynamics calculation"
 572          if (lang.eq.1) then
 573           write (iout,'(a/)')
 574      &      " with direct integration of Langevin equations"
 575          else if (lang.eq.2) then
 576           write (iout,'(a/)') " with TINKER stochasic MD integrator"
 577          else if (lang.eq.3) then
 578           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
 579          else if (lang.eq.4) then
 580           write (iout,'(a/)') " in overdamped mode"
 581          else
 582           write (iout,'(//a,i5)')
 583      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
 584           stop
 585          endif
 586          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 587          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
 588          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
 589          write (iout,'(a60,f10.5)')
 590      &   "Scaling factor of the friction forces:",scal_fric
 591          if (surfarea) write (iout,'(2a,i10,a)')
 592      &     "Friction coefficients will be scaled by solvent-accessible",
 593      &     " surface area every",reset_fricmat," steps."
 594         endif
 595 c Calculate friction coefficients and bounds of stochastic forces
 596         eta=6*pi*cPoise*etawat
 597         if(me.eq.king.or..not.out1file)
 598      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
 599      &   ,eta
 600         gamp=scal_fric*(pstok+rwat)*eta
 601         stdfp=dsqrt(2*Rb*t_bath/d_time)
 602         do i=1,ntyp
 603           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
 604           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 605         enddo
 606         if(me.eq.king.or..not.out1file)then
 607          write (iout,'(/2a/)')
 608      &   "Radii of site types and friction coefficients and std's of",
 609      &   " stochastic forces of fully exposed sites"
 610          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
 611          do i=1,ntyp
 612           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
 613      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
 614          enddo
 615         endif
 616       else if (tbf) then
 617         if(me.eq.king.or..not.out1file)then
 618          write (iout,'(a)') "Berendsen bath calculation"
 619          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 620          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
 621          if (reset_moment)
 622      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 623      &   count_reset_moment," steps"
 624          if (reset_vel)
 625      &    write (iout,'(a,i10,a)')
 626      &    "Velocities will be reset at random every",count_reset_vel,
 627      &   " steps"
 628         endif
 629       else
 630         if(me.eq.king.or..not.out1file)
 631      &   write (iout,'(a31)') "Microcanonical mode calculation"
 632       endif
 633       if(me.eq.king.or..not.out1file)then
 634        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
 635        if (usampl) then
 636           write(iout,*) "MD running with constraints."
 637           write(iout,*) "Equilibration time ", eq_time, " mtus."
 638           write(iout,*) "Constraining ", nfrag," fragments."
 639           write(iout,*) "Length of each fragment, weight and q0:"
 640           do iset=1,nset
 641            write (iout,*) "Set of restraints #",iset
 642            do i=1,nfrag
 643               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
 644      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
 645            enddo
 646            write(iout,*) "constraints between ", npair, "fragments."
 647            write(iout,*) "constraint pairs, weights and q0:"
 648            do i=1,npair
 649             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
 650      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
 651            enddo
 652            write(iout,*) "angle constraints within ", nfrag_back,
 653      &      "backbone fragments."
 654            if (loc_qlike) then
 655            write(iout,*) "fragment, weights, q0:"
 656            do i=1,nfrag_back
 657             write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
 658      &         ifrag_back(2,i,iset),
 659      &         wfrag_back(1,i,iset),qin_back(1,i,iset),
 660      &         wfrag_back(2,i,iset),qin_back(2,i,iset),
 661      &         wfrag_back(3,i,iset),qin_back(3,i,iset)
 662            enddo
 663            else
 664            write(iout,*) "fragment, weights:"
 665            do i=1,nfrag_back
 666             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
 667      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
 668      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
 669            enddo
 670            endif
 671           enddo
 672         iset=mod(kolor,nset)+1
 673        endif
 674       endif
 675       if(me.eq.king.or..not.out1file)
 676      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
 677       return
 678       end
 679 c------------------------------------------------------------------------------
 680       subroutine molread
 681 C
 682 C Read molecular data.
 683 C
 684       implicit real*8 (a-h,o-z)
 685       include 'DIMENSIONS'
 686 #ifdef MPI
 687       include 'mpif.h'
 688       integer error_msg
 689 #endif
 690       include 'COMMON.IOUNITS'
 691       include 'COMMON.GEO'
 692       include 'COMMON.VAR'
 693       include 'COMMON.INTERACT'
 694       include 'COMMON.LOCAL'
 695       include 'COMMON.NAMES'
 696       include 'COMMON.CHAIN'
 697       include 'COMMON.FFIELD'
 698       include 'COMMON.SBRIDGE'
 699       include 'COMMON.HEADER'
 700       include 'COMMON.CONTROL'
 701       include 'COMMON.DBASE'
 702       include 'COMMON.THREAD'
 703       include 'COMMON.CONTACTS'
 704       include 'COMMON.TORCNSTR'
 705       include 'COMMON.TIME1'
 706       include 'COMMON.BOUNDS'
 707       include 'COMMON.MD'
 708       include 'COMMON.SETUP'
 709       include 'COMMON.SHIELD'
 710       character*4 sequence(maxres)
 711       integer rescode
 712       double precision x(maxvar)
 713       character*256 pdbfile
 714       character*400 weightcard
 715       character*80 weightcard_t,ucase
 716       dimension itype_pdb(maxres)
 717       common /pizda/ itype_pdb
 718       logical seq_comp,fail
 719       double precision energia(0:n_ene)
 720       double precision secprob(3,maxdih_constr)
 721       integer ilen
 722       external ilen
 723       integer tperm
 724 C
 725 C Read PDB structure if applicable
 726 C
 727       if (indpdb.gt.0 .or. pdbref) then
 728         read(inp,'(a)') pdbfile
 729         if(me.eq.king.or..not.out1file)
 730      &   write (iout,'(2a)') 'PDB data will be read from file ',
 731      &   pdbfile(:ilen(pdbfile))
 732         open(ipdbin,file=pdbfile,status='old',err=33)
 733         goto 34
 734   33    write (iout,'(a)') 'Error opening PDB file.'
 735         stop
 736   34    continue
 737         call readpdb
 738         do i=1,2*nres
 739           do j=1,3
 740             crefjlee(j,i)=c(j,i)
 741           enddo
 742         enddo
 743 #ifdef DEBUG
 744         do i=1,nres
 745           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 746      &      (crefjlee(j,i+nres),j=1,3)
 747         enddo
 748 #endif
 749         if(me.eq.king.or..not.out1file)
 750      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 751      &   ' nstart_sup=',nstart_sup
 752         do i=1,nres
 753           itype_pdb(i)=itype(i)
 754         enddo
 755         close (ipdbin)
 756         nnt=nstart_sup
 757         nct=nstart_sup+nsup-1
 758         call contact(.false.,ncont_ref,icont_ref,co)
 759
 760         if (sideadd) then
 761          if(me.eq.king.or..not.out1file)
 762      &    write(iout,*)'Adding sidechains'
 763          maxsi=1000
 764          do i=2,nres-1
 765           iti=itype(i)
 766           if (iti.ne.10 .and. itype(i).ne.ntyp1) then
 767             nsi=0
 768             fail=.true.
 769             do while (fail.and.nsi.le.maxsi)
 770               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 771               nsi=nsi+1
 772             enddo
 773             if(fail) write(iout,*)'Adding sidechain failed for res ',
 774      &              i,' after ',nsi,' trials'
 775           endif
 776          enddo
 777         endif
 778       endif
 779       if (indpdb.eq.0) then
 780 C Read sequence if not taken from the pdb file.
 781         read (inp,*) nres
 782 c        print *,'nres=',nres
 783         if (iscode.gt.0) then
 784           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 785         else
 786           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 787         endif
 788 C Convert sequence to numeric code
 789         do i=1,nres
 790           itype(i)=rescode(i,sequence(i),iscode)
 791         enddo
 792 C Assign initial virtual bond lengths
 793 c        do i=2,nres
 794 c          vbld(i)=vbl
 795 c          vbld_inv(i)=vblinv
 796 c        enddo
 797 c        if (itype(1).eq.ntyp1) then
 798 c          vbld(2)=vbld(2)/2
 799 c          vbld_inv(2)=vbld_inv(2)*2
 800 c        endif
 801 c        if (itype(nres).eq.ntyp1) then
 802 c          vbld(nres)=vbld(nres)/2
 803 c          vbld_inv(nres)=vbld_inv(nres)*2
 804 c        endif
 805 c        do i=2,nres-1
 806 c          vbld(i+nres)=dsc(iabs(itype(i)))
 807 c          vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
 808 c          write (iout,*) "i",i," itype",itype(i),
 809 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 810 c        enddo
 811       endif
 812 c      print *,nres
 813 c      print '(20i4)',(itype(i),i=1,nres)
 814       do i=1,nres
 815 #ifdef PROCOR
 816         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 817 #else
 818         if (itype(i).eq.ntyp1) then
 819 #endif
 820           itel(i)=0
 821 #ifdef PROCOR
 822         else if (iabs(itype(i+1)).ne.20) then
 823 #else
 824         else if (iabs(itype(i)).ne.20) then
 825 #endif
 826           itel(i)=1
 827         else
 828           itel(i)=2
 829         endif
 830       enddo
 831       if(me.eq.king.or..not.out1file)then
 832        write (iout,*) "ITEL"
 833        do i=1,nres-1
 834          write (iout,*) i,itype(i),itel(i)
 835        enddo
 836        print *,'Call Read_Bridge.'
 837       endif
 838       nnt=1
 839       nct=nres
 840 cd      print *,'NNT=',NNT,' NCT=',NCT
 841       call seq2chains(nres,itype,nchain,chain_length,chain_border,
 842      &  ireschain)
 843       write(iout,*) "nres",nres," nchain",nchain
 844       do i=1,nchain
 845         write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
 846      &    chain_border(2,i)
 847       enddo
 848       call chain_symmetry(nchain,nres,itype,chain_border,
 849      &    chain_length,npermchain,tabpermchain)
 850 c      do i=1,nres
 851 c        write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
 852 c     &    ii=1,npermchain)
 853 c      enddo
 854       write(iout,*) "residue permutations"
 855       do i=1,nres
 856         write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
 857       enddo
 858       if (itype(1).eq.ntyp1) nnt=2
 859       if (itype(nres).eq.ntyp1) nct=nct-1
 860 #ifdef DFA
 861       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 862      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 863        call init_dfa_vars
 864        print*, 'init_dfa_vars finished!'
 865        call read_dfa_info
 866        print*, 'read_dfa_info finished!'
 867       endif
 868 #endif
 869       if (pdbref) then
 870         if(me.eq.king.or..not.out1file)
 871      &   write (iout,'(a,i3)') 'nsup=',nsup
 872         nstart_seq=nnt
 873         if (nsup.le.(nct-nnt+1)) then
 874           do i=0,nct-nnt+1-nsup
 875             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 876               nstart_seq=nnt+i
 877               goto 111
 878             endif
 879           enddo
 880           write (iout,'(a)')
 881      &            'Error - sequences to be superposed do not match.'
 882           stop
 883         else
 884           do i=0,nsup-(nct-nnt+1)
 885             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 886      &      then
 887               nstart_sup=nstart_sup+i
 888               nsup=nct-nnt+1
 889               goto 111
 890             endif
 891           enddo
 892           write (iout,'(a)')
 893      &            'Error - sequences to be superposed do not match.'
 894         endif
 895   111   continue
 896         if (nsup.eq.0) nsup=nct-nnt
 897         if (nstart_sup.eq.0) nstart_sup=nnt
 898         if (nstart_seq.eq.0) nstart_seq=nnt
 899         if(me.eq.king.or..not.out1file)
 900      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 901      &                 ' nstart_seq=',nstart_seq
 902       endif
 903       call read_bridge
 904 C 8/13/98 Set limits to generating the dihedral angles
 905       do i=1,nres
 906         phibound(1,i)=-pi
 907         phibound(2,i)=pi
 908       enddo
 909       read (inp,*) ndih_constr
 910       write (iout,*) "ndih_constr",ndih_constr
 911       if (ndih_constr.gt.0) then
 912         raw_psipred=.false.
 913 C        read (inp,*) ftors
 914         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
 915      &  i=1,ndih_constr)
 916 #ifdef MPI
 917         if(me.eq.king.or..not.out1file)then
 918 #endif
 919          write (iout,*)
 920      &   'There are',ndih_constr,' restraints on gamma angles.'
 921          do i=1,ndih_constr
 922           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
 923      &    ftors(i)
 924          enddo
 925 #ifdef MPI
 926         endif
 927 #endif
 928         do i=1,ndih_constr
 929           phi0(i)=deg2rad*phi0(i)
 930           drange(i)=deg2rad*drange(i)
 931         enddo
 932 C        if(me.eq.king.or..not.out1file)
 933 C     &   write (iout,*) 'FTORS',ftors
 934         do i=1,ndih_constr
 935           ii = idih_constr(i)
 936           phibound(1,ii) = phi0(i)-drange(i)
 937           phibound(2,ii) = phi0(i)+drange(i)
 938         enddo
 939 #ifdef MPI
 940         if (me.eq.king .or. .not.out1file) then
 941 #endif
 942          write (iout,'(a)') 'Boundaries in gamma angle sampling:'
 943          do i=1,nres
 944            write (iout,'(a3,i5,2f10.1)')
 945      &    restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 946          enddo
 947 #ifdef MP
 948         endif
 949 #endif
 950       else if (ndih_constr.lt.0) then
 951         raw_psipred=.true.
 952         call card_concat(weightcard)
 953         call reada(weightcard,"PHIHEL",phihel,50.0D0)
 954         call reada(weightcard,"PHIBET",phibet,180.0D0)
 955         call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
 956         call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
 957         call reada(weightcard,"WDIHC",wdihc,0.591D0)
 958         write (iout,*) "Weight of dihedral angle restraints",wdihc
 959         read(inp,'(9x,3f7.3)')
 960      &     (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
 961         write (iout,*) "The secprob array"
 962         do i=nnt,nct
 963           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
 964         enddo
 965         ndih_constr=0
 966         do i=nnt+3,nct
 967           if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
 968      &    .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
 969             ndih_constr=ndih_constr+1
 970             idih_constr(ndih_constr)=i
 971             sumv=0.0d0
 972             do j=1,3
 973               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
 974               sumv=sumv+vpsipred(j,ndih_constr)
 975             enddo
 976             do j=1,3
 977               vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
 978             enddo
 979             phibound(1,ndih_constr)=phihel*deg2rad
 980             phibound(2,ndih_constr)=phibet*deg2rad
 981             sdihed(1,ndih_constr)=sigmahel*deg2rad
 982             sdihed(2,ndih_constr)=sigmabet*deg2rad
 983           endif
 984         enddo
 985 #ifdef MPI
 986         if(me.eq.king.or..not.out1file)then
 987 #endif
 988          write (iout,*)
 989      &   'There are',ndih_constr,
 990      &   ' bimodal restraints on gamma angles.'
 991          do i=1,ndih_constr
 992            write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
 993      &      restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
 994      &      restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
 995      &      phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
 996      &      phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
 997      &      (vpsipred(j,i),j=1,3)
 998          enddo
 999 #ifdef MPI
1000         endif
1001 #endif
1002 c Raw psipred input
1003       endif
1004 C first setting the theta boundaries to 0 to pi
1005 C this mean that there is no energy penalty for any angle occuring this can be applied
1006 C for generate random conformation but is not implemented in this way
1007 C      do i=1,nres
1008 C        thetabound(1,i)=0
1009 C        thetabound(2,i)=pi
1010 C      enddo
1011 C begin reading theta constrains this is quartic constrains allowing to
1012 C have smooth second derivative
1013       if (with_theta_constr) then
1014 C with_theta_constr is keyword allowing for occurance of theta constrains
1015       read (inp,*) ntheta_constr
1016 C ntheta_constr is the number of theta constrains
1017       if (ntheta_constr.gt.0) then
1018 C        read (inp,*) ftors
1019         read (inp,*) (itheta_constr(i),theta_constr0(i),
1020      &  theta_drange(i),for_thet_constr(i),
1021      &  i=1,ntheta_constr)
1022 C the above code reads from 1 to ntheta_constr
1023 C itheta_constr(i) residue i for which is theta_constr
1024 C theta_constr0 the global minimum value
1025 C theta_drange is range for which there is no energy penalty
1026 C for_thet_constr is the force constant for quartic energy penalty
1027 C E=k*x**4
1028         if(me.eq.king.or..not.out1file)then
1029          write (iout,*)
1030      &   'There are',ntheta_constr,' constraints on phi angles.'
1031          do i=1,ntheta_constr
1032           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1033      &    theta_drange(i),
1034      &    for_thet_constr(i)
1035          enddo
1036         endif
1037         do i=1,ntheta_constr
1038           theta_constr0(i)=deg2rad*theta_constr0(i)
1039           theta_drange(i)=deg2rad*theta_drange(i)
1040         enddo
1041 C        if(me.eq.king.or..not.out1file)
1042 C     &   write (iout,*) 'FTORS',ftors
1043 C        do i=1,ntheta_constr
1044 C          ii = itheta_constr(i)
1045 C          thetabound(1,ii) = phi0(i)-drange(i)
1046 C          thetabound(2,ii) = phi0(i)+drange(i)
1047 C        enddo
1048       endif ! ntheta_constr.gt.0
1049       endif! with_theta_constr
1050 C
1051 C      with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1052 C      write (iout,*) "with_dihed_constr ",with_dihed_constr
1053 c--- Zscore rms -------
1054       if (nz_start.eq.0) nz_start=nnt
1055       if (nz_end.eq.0 .and. nsup.gt.0) then
1056         nz_end=nnt+nsup-1
1057       else if (nz_end.eq.0) then
1058         nz_end=nct
1059       endif
1060       if(me.eq.king.or..not.out1file)then
1061        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1062        write (iout,*) 'IZ_SC=',iz_sc
1063       endif
1064 c----------------------
1065       call init_int_table
1066       if (refstr) then
1067         if (.not.pdbref) then
1068           call read_angles(inp,*38)
1069           call bond_regular
1070           goto 39
1071    38     write (iout,'(a)') 'Error reading reference structure.'
1072 #ifdef MPI
1073           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1074           stop 'Error reading reference structure'
1075 #endif
1076    39     call chainbuild_extconf
1077           call setup_var
1078 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1079           nstart_sup=nnt
1080           nstart_seq=nnt
1081           nsup=nct-nnt+1
1082           do i=1,2*nres
1083             do j=1,3
1084               cref(j,i)=c(j,i)
1085             enddo
1086           enddo
1087           call contact(.true.,ncont_ref,icont_ref,co)
1088         endif
1089       endif
1090 c        print *, "A TU"
1091       write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1092       call flush(iout)
1093       if (constr_dist.gt.0) call read_dist_constr
1094       write (iout,*) "After read_dist_constr nhpb",nhpb
1095       if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1096       call hpb_partition
1097       call NMRpeak_partition
1098       if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1099       if (pdbref) then
1100         if(me.eq.king.or..not.out1file)
1101      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1102         do i=1,ncont_ref
1103           do j=1,2
1104             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1105           enddo
1106           if(me.eq.king.or..not.out1file)
1107      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1108      &     icont_ref(1,i),' ',
1109      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1110         enddo
1111       endif
1112       write (iout,*) "calling read_saxs_consrtr",nsaxs
1113       if (nsaxs.gt.0) call read_saxs_constr
1114
1115       if (constr_homology.gt.0) then
1116         call read_constr_homology
1117         if (indpdb.gt.0 .or. pdbref) then
1118           do i=1,2*nres
1119             do j=1,3
1120               c(j,i)=crefjlee(j,i)
1121               cref(j,i)=crefjlee(j,i)
1122             enddo
1123           enddo
1124         endif
1125 #ifdef DEBUG
1126         write (iout,*) "sc_loc_geom: Array C"
1127         do i=1,nres
1128           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1129      &      (c(j,i+nres),j=1,3)
1130         enddo
1131         write (iout,*) "Array Cref"
1132         do i=1,nres
1133           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1134      &      (cref(j,i+nres),j=1,3)
1135         enddo
1136 #endif
1137        call int_from_cart1(.false.)
1138        call sc_loc_geom(.false.)
1139        do i=1,nres
1140          thetaref(i)=theta(i)
1141          phiref(i)=phi(i)
1142        enddo
1143        do i=1,nres-1
1144          do j=1,3
1145            dc(j,i)=c(j,i+1)-c(j,i)
1146            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1147          enddo
1148        enddo
1149        do i=2,nres-1
1150          do j=1,3
1151            dc(j,i+nres)=c(j,i+nres)-c(j,i)
1152            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1153          enddo
1154        enddo
1155       else
1156         homol_nset=0
1157       endif
1158
1159
1160 C      endif
1161       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1162      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1163      &    modecalc.ne.10) then
1164 C If input structure hasn't been supplied from the PDB file read or generate
1165 C initial geometry.
1166         if (iranconf.eq.0 .and. .not. extconf) then
1167           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1168      &     write (iout,'(a)') 'Initial geometry will be read in.'
1169           if (read_cart) then
1170             read(inp,'(8f10.5)',end=36,err=36)
1171      &       ((c(l,k),l=1,3),k=1,nres),
1172      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1173             write (iout,*) "Exit READ_CART"
1174 c            write (iout,'(8f10.5)')
1175 c     &       ((c(l,k),l=1,3),k=1,nres),
1176 c     &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1177             call cartprint
1178             do i=1,nres-1
1179               do j=1,3
1180                 dc(j,i)=c(j,i+1)-c(j,i)
1181 c                dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1182               enddo
1183             enddo
1184             do i=nnt,nct
1185               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1186                 do j=1,3
1187                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1188 c                  dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1189                 enddo
1190               else
1191                 do j=1,3
1192                   dc(j,i+nres)=0.0d0
1193 c                  dc_norm(j,i+nres)=0.0d0
1194                 enddo
1195               endif
1196             enddo
1197             call int_from_cart1(.true.)
1198             write (iout,*) "Finish INT_TO_CART"
1199 c      write (iout,*) "DC"
1200 c      do i=1,nres
1201 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1202 c     &    (dc(j,i+nres),j=1,3)
1203 c      enddo
1204 c            call cartprint
1205 c            call setup_var
1206 c            return
1207           else
1208             call read_angles(inp,*36)
1209             call bond_regular
1210             call chainbuild_extconf
1211           endif
1212           goto 37
1213    36     write (iout,'(a)') 'Error reading angle file.'
1214 #ifdef MPI
1215           call mpi_finalize( MPI_COMM_WORLD,IERR )
1216 #endif
1217           stop 'Error reading angle file.'
1218    37     continue
1219         else if (extconf) then
1220           if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1221      &      write (iout,'(a)') 'Extended chain initial geometry.'
1222           do i=3,nres
1223             theta(i)=90d0*deg2rad
1224           enddo
1225           do i=4,nres
1226             phi(i)=180d0*deg2rad
1227           enddo
1228           do i=2,nres-1
1229             alph(i)=110d0*deg2rad
1230           enddo
1231           do i=2,nres-1
1232             omeg(i)=-120d0*deg2rad
1233             if (itype(i).le.0) omeg(i)=-omeg(i)
1234           enddo
1235           call bond_regular
1236           call chainbuild_extconf
1237         else
1238           if(me.eq.king.or..not.out1file)
1239      &     write (iout,'(a)') 'Random-generated initial geometry.'
1240           call bond_regular
1241 #ifdef MPI
1242           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1243      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1244 #endif
1245             do itrial=1,100
1246               itmp=1
1247               call gen_rand_conf(itmp,*30)
1248               goto 40
1249    30         write (iout,*) 'Failed to generate random conformation',
1250      &          ', itrial=',itrial
1251               write (*,*) 'Processor:',me,
1252      &          ' Failed to generate random conformation',
1253      &          ' itrial=',itrial
1254               call intout
1255
1256 #ifdef AIX
1257               call flush_(iout)
1258 #else
1259               call flush(iout)
1260 #endif
1261             enddo
1262             write (iout,'(a,i3,a)') 'Processor:',me,
1263      &        ' error in generating random conformation.'
1264             write (*,'(a,i3,a)') 'Processor:',me,
1265      &        ' error in generating random conformation.'
1266             call flush(iout)
1267 #ifdef MPI
1268             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1269    40       continue
1270           endif
1271 #else
1272           write (*,'(a)')
1273      &      ' error in generating random conformation.'
1274           stop
1275    40     continue
1276 #endif
1277         endif
1278       elseif (modecalc.eq.4) then
1279         read (inp,'(a)') intinname
1280         open (intin,file=intinname,status='old',err=333)
1281         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1282      &  write (iout,'(a)') 'intinname',intinname
1283         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1284         goto 334
1285   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1286 #ifdef MPI
1287         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1288 #endif
1289         stop 'Error opening angle file.'
1290   334   continue
1291
1292       endif
1293 C Generate distance constraints, if the PDB structure is to be regularized.
1294       if (nthread.gt.0) then
1295         call read_threadbase
1296       endif
1297       call setup_var
1298       if (me.eq.king .or. .not. out1file)
1299      & call intout
1300       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1301         write (iout,'(/a,i3,a)')
1302      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1303         write (iout,'(20i4)') (iss(i),i=1,ns)
1304        if (dyn_ss) then
1305           write(iout,*)"Running with dynamic disulfide-bond formation"
1306        else
1307         write (iout,'(/a/)') 'Pre-formed links are:'
1308         do i=1,nss
1309           i1=ihpb(i)-nres
1310           i2=jhpb(i)-nres
1311           it1=itype(i1)
1312           it2=itype(i2)
1313           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1314      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1315      &    ebr,forcon(i)
1316         enddo
1317         write (iout,'(a)')
1318        endif
1319       endif
1320       if (ns.gt.0.and.dyn_ss) then
1321           do i=nss+1,nhpb
1322             ihpb(i-nss)=ihpb(i)
1323             jhpb(i-nss)=jhpb(i)
1324             forcon(i-nss)=forcon(i)
1325             dhpb(i-nss)=dhpb(i)
1326           enddo
1327           nhpb=nhpb-nss
1328           nss=0
1329           call hpb_partition
1330           do i=1,ns
1331             dyn_ss_mask(iss(i))=.true.
1332           enddo
1333       endif
1334 c --- test
1335 c      call cartprint
1336 c      write (iout,*) "DC"
1337 c      do i=1,nres
1338 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1339 c     &    (dc(j,i+nres),j=1,3)
1340 c      enddo
1341       if (i2ndstr.gt.0) call secstrp2dihc
1342 c      call geom_to_var(nvar,x)
1343 c      call etotal(energia(0))
1344 c      call enerprint(energia(0))
1345 c      call briefout(0,etot)
1346 c      stop
1347 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1348 cd    write (iout,'(a)') 'Variable list:'
1349 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1350 #ifdef MPI
1351       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1352      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1353      &  'Processor',myrank,': end reading molecular data.'
1354 #endif
1355 c      print *,"A TU?"
1356       return
1357       end
1358 c--------------------------------------------------------------------------
1359       logical function seq_comp(itypea,itypeb,length)
1360       implicit none
1361       integer length,itypea(length),itypeb(length)
1362       integer i
1363       do i=1,length
1364         if (itypea(i).ne.itypeb(i)) then
1365           seq_comp=.false.
1366           return
1367         endif
1368       enddo
1369       seq_comp=.true.
1370       return
1371       end
1372 c-----------------------------------------------------------------------------
1373       subroutine read_bridge
1374 C Read information about disulfide bridges.
1375       implicit real*8 (a-h,o-z)
1376       include 'DIMENSIONS'
1377 #ifdef MPI
1378       include 'mpif.h'
1379 #endif
1380       include 'COMMON.IOUNITS'
1381       include 'COMMON.GEO'
1382       include 'COMMON.VAR'
1383       include 'COMMON.INTERACT'
1384       include 'COMMON.LOCAL'
1385       include 'COMMON.NAMES'
1386       include 'COMMON.CHAIN'
1387       include 'COMMON.FFIELD'
1388       include 'COMMON.SBRIDGE'
1389       include 'COMMON.HEADER'
1390       include 'COMMON.CONTROL'
1391       include 'COMMON.DBASE'
1392       include 'COMMON.THREAD'
1393       include 'COMMON.TIME1'
1394       include 'COMMON.SETUP'
1395 C Read bridging residues.
1396       read (inp,*) ns,(iss(i),i=1,ns)
1397       print *,'ns=',ns
1398       if(me.eq.king.or..not.out1file)
1399      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1400 C Check whether the specified bridging residues are cystines.
1401       do i=1,ns
1402         if (itype(iss(i)).ne.1) then
1403           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1404      &   'Do you REALLY think that the residue ',
1405      &    restyp(itype(iss(i))),i,
1406      &   ' can form a disulfide bridge?!!!'
1407           write (*,'(2a,i3,a)')
1408      &   'Do you REALLY think that the residue ',
1409      &    restyp(itype(iss(i))),i,
1410      &   ' can form a disulfide bridge?!!!'
1411 #ifdef MPI
1412          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1413          stop
1414 #endif
1415         endif
1416       enddo
1417 C Read preformed bridges.
1418       if (ns.gt.0) then
1419       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1420       if(fg_rank.eq.0)
1421      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1422       if (nss.gt.0) then
1423         nhpb=nss
1424 C Check if the residues involved in bridges are in the specified list of
1425 C bridging residues.
1426         do i=1,nss
1427           do j=1,i-1
1428             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1429      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1430               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1431      &      ' contains residues present in other pairs.'
1432               write (*,'(a,i3,a)') 'Disulfide pair',i,
1433      &      ' contains residues present in other pairs.'
1434 #ifdef MPI
1435               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1436               stop
1437 #endif
1438             endif
1439           enddo
1440           do j=1,ns
1441             if (ihpb(i).eq.iss(j)) goto 10
1442           enddo
1443           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1444    10     continue
1445           do j=1,ns
1446             if (jhpb(i).eq.iss(j)) goto 20
1447           enddo
1448           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1449    20     continue
1450           dhpb(i)=dbr
1451           forcon(i)=fbr
1452         enddo
1453         do i=1,nss
1454           ihpb(i)=ihpb(i)+nres
1455           jhpb(i)=jhpb(i)+nres
1456         enddo
1457       endif
1458       endif
1459       return
1460       end
1461 c----------------------------------------------------------------------------
1462       subroutine read_x(kanal,*)
1463       implicit real*8 (a-h,o-z)
1464       include 'DIMENSIONS'
1465       include 'COMMON.GEO'
1466       include 'COMMON.VAR'
1467       include 'COMMON.CHAIN'
1468       include 'COMMON.IOUNITS'
1469       include 'COMMON.CONTROL'
1470       include 'COMMON.LOCAL'
1471       include 'COMMON.INTERACT'
1472 c Read coordinates from input
1473 c
1474       read(kanal,'(8f10.5)',end=10,err=10)
1475      &  ((c(l,k),l=1,3),k=1,nres),
1476      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1477       do j=1,3
1478         c(j,nres+1)=c(j,1)
1479         c(j,2*nres)=c(j,nres)
1480       enddo
1481       call int_from_cart1(.false.)
1482       do i=1,nres-1
1483         do j=1,3
1484           dc(j,i)=c(j,i+1)-c(j,i)
1485           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1486         enddo
1487       enddo
1488       do i=nnt,nct
1489         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1490           do j=1,3
1491             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1492             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1493           enddo
1494         endif
1495       enddo
1496
1497       return
1498    10 return1
1499       end
1500 c----------------------------------------------------------------------------
1501       subroutine read_threadbase
1502       implicit real*8 (a-h,o-z)
1503       include 'DIMENSIONS'
1504       include 'COMMON.IOUNITS'
1505       include 'COMMON.GEO'
1506       include 'COMMON.VAR'
1507       include 'COMMON.INTERACT'
1508       include 'COMMON.LOCAL'
1509       include 'COMMON.NAMES'
1510       include 'COMMON.CHAIN'
1511       include 'COMMON.FFIELD'
1512       include 'COMMON.SBRIDGE'
1513       include 'COMMON.HEADER'
1514       include 'COMMON.CONTROL'
1515       include 'COMMON.DBASE'
1516       include 'COMMON.THREAD'
1517       include 'COMMON.TIME1'
1518 C Read pattern database for threading.
1519       read (icbase,*) nseq
1520       do i=1,nseq
1521         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1522      &   nres_base(2,i),nres_base(3,i)
1523         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1524      &   nres_base(1,i))
1525 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1526 c    &   nres_base(2,i),nres_base(3,i)
1527 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1528 c    &   nres_base(1,i))
1529       enddo
1530       close (icbase)
1531       if (weidis.eq.0.0D0) weidis=0.1D0
1532       do i=nnt,nct
1533         do j=i+2,nct
1534           nhpb=nhpb+1
1535           ihpb(nhpb)=i
1536           jhpb(nhpb)=j
1537           forcon(nhpb)=weidis
1538         enddo
1539       enddo
1540       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1541       write (iout,'(a,i5)') 'nexcl: ',nexcl
1542       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1543       return
1544       end
1545 c------------------------------------------------------------------------------
1546       subroutine setup_var
1547       implicit real*8 (a-h,o-z)
1548       include 'DIMENSIONS'
1549       include 'COMMON.IOUNITS'
1550       include 'COMMON.GEO'
1551       include 'COMMON.VAR'
1552       include 'COMMON.INTERACT'
1553       include 'COMMON.LOCAL'
1554       include 'COMMON.NAMES'
1555       include 'COMMON.CHAIN'
1556       include 'COMMON.FFIELD'
1557       include 'COMMON.SBRIDGE'
1558       include 'COMMON.HEADER'
1559       include 'COMMON.CONTROL'
1560       include 'COMMON.DBASE'
1561       include 'COMMON.THREAD'
1562       include 'COMMON.TIME1'
1563 C Set up variable list.
1564       ntheta=nres-2
1565       nphi=nres-3
1566       nvar=ntheta+nphi
1567       nside=0
1568       write (iout,*) "SETUP_VAR ialph"
1569       do i=2,nres-1
1570         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1571           nside=nside+1
1572           ialph(i,1)=nvar+nside
1573           ialph(nside,2)=i
1574         endif
1575       enddo
1576       if (indphi.gt.0) then
1577         nvar=nphi
1578       else if (indback.gt.0) then
1579         nvar=nphi+ntheta
1580       else
1581         nvar=nvar+2*nside
1582       endif
1583       write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1584       return
1585       end
1586 c----------------------------------------------------------------------------
1587       subroutine gen_dist_constr
1588 C Generate CA distance constraints.
1589       implicit real*8 (a-h,o-z)
1590       include 'DIMENSIONS'
1591       include 'COMMON.IOUNITS'
1592       include 'COMMON.GEO'
1593       include 'COMMON.VAR'
1594       include 'COMMON.INTERACT'
1595       include 'COMMON.LOCAL'
1596       include 'COMMON.NAMES'
1597       include 'COMMON.CHAIN'
1598       include 'COMMON.FFIELD'
1599       include 'COMMON.SBRIDGE'
1600       include 'COMMON.HEADER'
1601       include 'COMMON.CONTROL'
1602       include 'COMMON.DBASE'
1603       include 'COMMON.THREAD'
1604       include 'COMMON.TIME1'
1605       dimension itype_pdb(maxres)
1606       common /pizda/ itype_pdb
1607       character*2 iden
1608 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1609 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1610 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1611 cd     & ' nsup',nsup
1612       do i=nstart_sup,nstart_sup+nsup-1
1613 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1614 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1615         do j=i+2,nstart_sup+nsup-1
1616           nhpb=nhpb+1
1617           ihpb(nhpb)=i+nstart_seq-nstart_sup
1618           jhpb(nhpb)=j+nstart_seq-nstart_sup
1619           forcon(nhpb)=weidis
1620           dhpb(nhpb)=dist(i,j)
1621         enddo
1622       enddo
1623 cd      write (iout,'(a)') 'Distance constraints:'
1624 cd      do i=nss+1,nhpb
1625 cd        ii=ihpb(i)
1626 cd        jj=jhpb(i)
1627 cd        iden='CA'
1628 cd        if (ii.gt.nres) then
1629 cd          iden='SC'
1630 cd          ii=ii-nres
1631 cd          jj=jj-nres
1632 cd        endif
1633 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1634 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1635 cd     &  dhpb(i),forcon(i)
1636 cd      enddo
1637       return
1638       end
1639 c----------------------------------------------------------------------------
1640       subroutine map_read
1641       implicit real*8 (a-h,o-z)
1642       include 'DIMENSIONS'
1643       include 'COMMON.MAP'
1644       include 'COMMON.IOUNITS'
1645       character*3 angid(4) /'THE','PHI','ALP','OME'/
1646       character*80 mapcard,ucase
1647       do imap=1,nmap
1648         read (inp,'(a)') mapcard
1649         mapcard=ucase(mapcard)
1650         if (index(mapcard,'PHI').gt.0) then
1651           kang(imap)=1
1652         else if (index(mapcard,'THE').gt.0) then
1653           kang(imap)=2
1654         else if (index(mapcard,'ALP').gt.0) then
1655           kang(imap)=3
1656         else if (index(mapcard,'OME').gt.0) then
1657           kang(imap)=4
1658         else
1659           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1660           stop 'Error - illegal variable spec in MAP card.'
1661         endif
1662         call readi (mapcard,'RES1',res1(imap),0)
1663         call readi (mapcard,'RES2',res2(imap),0)
1664         if (res1(imap).eq.0) then
1665           res1(imap)=res2(imap)
1666         else if (res2(imap).eq.0) then
1667           res2(imap)=res1(imap)
1668         endif
1669         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1670           write (iout,'(a)')
1671      &    'Error - illegal definition of variable group in MAP.'
1672           stop 'Error - illegal definition of variable group in MAP.'
1673         endif
1674         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1675         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1676         call readi(mapcard,'NSTEP',nstep(imap),0)
1677         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1678           write (iout,'(a)')
1679      &     'Illegal boundary and/or step size specification in MAP.'
1680           stop 'Illegal boundary and/or step size specification in MAP.'
1681         endif
1682       enddo ! imap
1683       return
1684       end
1685 c----------------------------------------------------------------------------
1686       subroutine csaread
1687       implicit real*8 (a-h,o-z)
1688       include 'DIMENSIONS'
1689       include 'COMMON.IOUNITS'
1690       include 'COMMON.GEO'
1691       include 'COMMON.CSA'
1692       include 'COMMON.BANK'
1693       include 'COMMON.CONTROL'
1694       character*80 ucase
1695       character*620 mcmcard
1696       call card_concat(mcmcard)
1697
1698       call readi(mcmcard,'NCONF',nconf,50)
1699       call readi(mcmcard,'NADD',nadd,0)
1700       call readi(mcmcard,'JSTART',jstart,1)
1701       call readi(mcmcard,'JEND',jend,1)
1702       call readi(mcmcard,'NSTMAX',nstmax,500000)
1703       call readi(mcmcard,'N0',n0,1)
1704       call readi(mcmcard,'N1',n1,6)
1705       call readi(mcmcard,'N2',n2,4)
1706       call readi(mcmcard,'N3',n3,0)
1707       call readi(mcmcard,'N4',n4,0)
1708       call readi(mcmcard,'N5',n5,0)
1709       call readi(mcmcard,'N6',n6,10)
1710       call readi(mcmcard,'N7',n7,0)
1711       call readi(mcmcard,'N8',n8,0)
1712       call readi(mcmcard,'N9',n9,0)
1713       call readi(mcmcard,'N14',n14,0)
1714       call readi(mcmcard,'N15',n15,0)
1715       call readi(mcmcard,'N16',n16,0)
1716       call readi(mcmcard,'N17',n17,0)
1717       call readi(mcmcard,'N18',n18,0)
1718
1719       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1720
1721       call readi(mcmcard,'NDIFF',ndiff,2)
1722       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1723       call readi(mcmcard,'IS1',is1,1)
1724       call readi(mcmcard,'IS2',is2,8)
1725       call readi(mcmcard,'NRAN0',nran0,4)
1726       call readi(mcmcard,'NRAN1',nran1,2)
1727       call readi(mcmcard,'IRR',irr,1)
1728       call readi(mcmcard,'NSEED',nseed,20)
1729       call readi(mcmcard,'NTOTAL',ntotal,10000)
1730       call reada(mcmcard,'CUT1',cut1,2.0d0)
1731       call reada(mcmcard,'CUT2',cut2,5.0d0)
1732       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1733       call readi(mcmcard,'ICMAX',icmax,3)
1734       call readi(mcmcard,'IRESTART',irestart,0)
1735 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1736       ntbankm=0
1737 c!bankt
1738       call reada(mcmcard,'DELE',dele,20.0d0)
1739       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1740       call readi(mcmcard,'IREF',iref,0)
1741       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1742       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1743       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1744       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1745       write (iout,*) "NCONF_IN",nconf_in
1746       return
1747       end
1748 c----------------------------------------------------------------------------
1749 cfmc      subroutine mcmfread
1750 cfmc      implicit real*8 (a-h,o-z)
1751 cfmc      include 'DIMENSIONS'
1752 cfmc      include 'COMMON.MCMF'
1753 cfmc      include 'COMMON.IOUNITS'
1754 cfmc      include 'COMMON.GEO'
1755 cfmc      character*80 ucase
1756 cfmc      character*620 mcmcard
1757 cfmc      call card_concat(mcmcard)
1758 cfmc
1759 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1760 cfmc      write(iout,*)'MAXRANT=',maxrant
1761 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1762 cfmc      write(iout,*)'MAXFAM=',maxfam
1763 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1764 cfmc      write(iout,*)'NNET1=',nnet1
1765 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1766 cfmc      write(iout,*)'NNET2=',nnet2
1767 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1768 cfmc      write(iout,*)'NNET3=',nnet3
1769 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1770 cfmc      write(iout,*)'ILASTT=',ilastt
1771 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1772 cfmc      write(iout,*)'MAXSTR=',maxstr
1773 cfmc      maxstr_f=maxstr/maxfam
1774 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1775 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1776 cfmc      write(iout,*)'NMCMF=',nmcmf
1777 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1778 cfmc      write(iout,*)'IFOCUS=',ifocus
1779 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1780 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1781 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1782 cfmc      write(iout,*)'INTPRT=',intprt
1783 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1784 cfmc      write(iout,*)'IPRT=',iprt
1785 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1786 cfmc      write(iout,*)'IMAXTR=',imaxtr
1787 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1788 cfmc      write(iout,*)'MAXEVEN=',maxeven
1789 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1790 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1791 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1792 cfmc      write(iout,*)'INIMIN=',inimin
1793 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1794 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1795 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1796 cfmc      write(iout,*)'NTHREAD=',nthread
1797 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1798 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1799 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1800 cfmc      write(iout,*)'MAXPERT=',maxpert
1801 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1802 cfmc      write(iout,*)'IRMSD=',irmsd
1803 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1804 cfmc      write(iout,*)'DENEMIN=',denemin
1805 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1806 cfmc      write(iout,*)'RCUT1S=',rcut1s
1807 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1808 cfmc      write(iout,*)'RCUT1E=',rcut1e
1809 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1810 cfmc      write(iout,*)'RCUT2S=',rcut2s
1811 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1812 cfmc      write(iout,*)'RCUT2E=',rcut2e
1813 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1814 cfmc      write(iout,*)'DPERT1=',d_pert1
1815 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1816 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1817 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1818 cfmc      write(iout,*)'DPERT2=',d_pert2
1819 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1820 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1821 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1822 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1823 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1824 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1825 cfmc      d_pert1=deg2rad*d_pert1
1826 cfmc      d_pert1a=deg2rad*d_pert1a
1827 cfmc      d_pert2=deg2rad*d_pert2
1828 cfmc      d_pert2a=deg2rad*d_pert2a
1829 cfmc      d_pert2b=deg2rad*d_pert2b
1830 cfmc      d_pert2c=deg2rad*d_pert2c
1831 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1832 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1833 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1834 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1835 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1836 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1837 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1838 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1839 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1840 cfmc      write(iout,*)'RCUTINI=',rcutini
1841 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1842 cfmc      write(iout,*)'GRAT=',grat
1843 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1844 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1845 cfmc
1846 cfmc      return
1847 cfmc      end
1848 c----------------------------------------------------------------------------
1849       subroutine mcmread
1850       implicit real*8 (a-h,o-z)
1851       include 'DIMENSIONS'
1852       include 'COMMON.MCM'
1853       include 'COMMON.MCE'
1854       include 'COMMON.IOUNITS'
1855       character*80 ucase
1856       character*320 mcmcard
1857       call card_concat(mcmcard)
1858       call readi(mcmcard,'MAXACC',maxacc,100)
1859       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1860       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1861       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1862       call readi(mcmcard,'MAXREPM',maxrepm,200)
1863       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1864       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1865       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1866       call reada(mcmcard,'E_UP',e_up,5.0D0)
1867       call reada(mcmcard,'DELTE',delte,0.1D0)
1868       call readi(mcmcard,'NSWEEP',nsweep,5)
1869       call readi(mcmcard,'NSTEPH',nsteph,0)
1870       call readi(mcmcard,'NSTEPC',nstepc,0)
1871       call reada(mcmcard,'TMIN',tmin,298.0D0)
1872       call reada(mcmcard,'TMAX',tmax,298.0D0)
1873       call readi(mcmcard,'NWINDOW',nwindow,0)
1874       call readi(mcmcard,'PRINT_MC',print_mc,0)
1875       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1876       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1877       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1878       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1879       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1880       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1881       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1882       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1883       if (nwindow.gt.0) then
1884         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1885         do i=1,nwindow
1886           winlen(i)=winend(i)-winstart(i)+1
1887         enddo
1888       endif
1889       if (tmax.lt.tmin) tmax=tmin
1890       if (tmax.eq.tmin) then
1891         nstepc=0
1892         nsteph=0
1893       endif
1894       if (nstepc.gt.0 .and. nsteph.gt.0) then
1895         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1896         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1897       endif
1898 C Probabilities of different move types
1899       sumpro_type(0)=0.0D0
1900       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1901       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1902       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1903       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1904       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1905       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1906       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1907       do i=1,MaxMoveType
1908         print *,'i',i,' sumprotype',sumpro_type(i)
1909         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1910         print *,'i',i,' sumprotype',sumpro_type(i)
1911       enddo
1912       return
1913       end
1914 c----------------------------------------------------------------------------
1915       subroutine read_minim
1916       implicit real*8 (a-h,o-z)
1917       include 'DIMENSIONS'
1918       include 'COMMON.MINIM'
1919       include 'COMMON.IOUNITS'
1920       character*80 ucase
1921       character*320 minimcard
1922       call card_concat(minimcard)
1923       call readi(minimcard,'MAXMIN',maxmin,2000)
1924       call readi(minimcard,'MAXFUN',maxfun,5000)
1925       call readi(minimcard,'MINMIN',minmin,maxmin)
1926       call readi(minimcard,'MINFUN',minfun,maxmin)
1927       call reada(minimcard,'TOLF',tolf,1.0D-2)
1928       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1929       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1930       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1931       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1932       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1933      &         'Options in energy minimization:'
1934       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1935      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1936      & 'MinMin:',MinMin,' MinFun:',MinFun,
1937      & ' TolF:',TolF,' RTolF:',RTolF
1938       return
1939       end
1940 c----------------------------------------------------------------------------
1941       subroutine read_angles(kanal,*)
1942       implicit real*8 (a-h,o-z)
1943       include 'DIMENSIONS'
1944       include 'COMMON.GEO'
1945       include 'COMMON.VAR'
1946       include 'COMMON.CHAIN'
1947       include 'COMMON.IOUNITS'
1948       include 'COMMON.CONTROL'
1949 c Read angles from input
1950 c
1951        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1952        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1953        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1954        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1955
1956        do i=1,nres
1957 c 9/7/01 avoid 180 deg valence angle
1958         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1959 c
1960         theta(i)=deg2rad*theta(i)
1961         phi(i)=deg2rad*phi(i)
1962         alph(i)=deg2rad*alph(i)
1963         omeg(i)=deg2rad*omeg(i)
1964        enddo
1965       return
1966    10 return1
1967       end
1968 c----------------------------------------------------------------------------
1969       subroutine reada(rekord,lancuch,wartosc,default)
1970       implicit none
1971       character*(*) rekord,lancuch
1972       double precision wartosc,default
1973       integer ilen,iread
1974       external ilen
1975       iread=index(rekord,lancuch)
1976       if (iread.eq.0) then
1977         wartosc=default
1978         return
1979       endif
1980       iread=iread+ilen(lancuch)+1
1981       read (rekord(iread:),*,err=10,end=10) wartosc
1982       return
1983   10  wartosc=default
1984       return
1985       end
1986 c----------------------------------------------------------------------------
1987       subroutine readi(rekord,lancuch,wartosc,default)
1988       implicit none
1989       character*(*) rekord,lancuch
1990       integer wartosc,default
1991       integer ilen,iread
1992       external ilen
1993       iread=index(rekord,lancuch)
1994       if (iread.eq.0) then
1995         wartosc=default
1996         return
1997       endif
1998       iread=iread+ilen(lancuch)+1
1999       read (rekord(iread:),*,err=10,end=10) wartosc
2000       return
2001   10  wartosc=default
2002       return
2003       end
2004 c----------------------------------------------------------------------------
2005       subroutine multreadi(rekord,lancuch,tablica,dim,default)
2006       implicit none
2007       integer dim,i
2008       integer tablica(dim),default
2009       character*(*) rekord,lancuch
2010       character*80 aux
2011       integer ilen,iread
2012       external ilen
2013       do i=1,dim
2014         tablica(i)=default
2015       enddo
2016       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2017       if (iread.eq.0) return
2018       iread=iread+ilen(lancuch)+1
2019       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2020    10 return
2021       end
2022 c----------------------------------------------------------------------------
2023       subroutine multreada(rekord,lancuch,tablica,dim,default)
2024       implicit none
2025       integer dim,i
2026       double precision tablica(dim),default
2027       character*(*) rekord,lancuch
2028       character*80 aux
2029       integer ilen,iread
2030       external ilen
2031       do i=1,dim
2032         tablica(i)=default
2033       enddo
2034       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2035       if (iread.eq.0) return
2036       iread=iread+ilen(lancuch)+1
2037       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2038    10 return
2039       end
2040 c----------------------------------------------------------------------------
2041       subroutine openunits
2042       implicit real*8 (a-h,o-z)
2043       include 'DIMENSIONS'
2044 #ifdef MPI
2045       include 'mpif.h'
2046       character*16 form,nodename
2047       integer nodelen
2048 #endif
2049       include 'COMMON.SETUP'
2050       include 'COMMON.IOUNITS'
2051       include 'COMMON.MD'
2052       include 'COMMON.CONTROL'
2053       integer lenpre,lenpot,ilen,lentmp
2054       external ilen
2055       character*3 out1file_text,ucase
2056       character*3 ll
2057       external ucase
2058       tmpdir=""
2059       out1file_text="YES"
2060 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2061       call getenv_loc("PREFIX",prefix)
2062       pref_orig = prefix
2063       call getenv_loc("POT",pot)
2064       call getenv_loc("DIRTMP",tmpdir)
2065       call getenv_loc("CURDIR",curdir)
2066       call getenv_loc("OUT1FILE",out1file_text)
2067 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2068       out1file_text=ucase(out1file_text)
2069       if (out1file_text(1:1).eq."Y") then
2070         out1file=.true.
2071       else
2072         out1file=fg_rank.gt.0
2073       endif
2074       lenpre=ilen(prefix)
2075       lenpot=ilen(pot)
2076       lentmp=ilen(tmpdir)
2077       if (lentmp.gt.0) then
2078           write (*,'(80(1h!))')
2079           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
2080           write (*,'(80(1h!))')
2081           write (*,*)"All output files will be on node /tmp directory."
2082 #ifdef MPI
2083         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2084         if (me.eq.king) then
2085           write (*,*) "The master node is ",nodename
2086         else if (fg_rank.eq.0) then
2087           write (*,*) "I am the CG slave node ",nodename
2088         else
2089           write (*,*) "I am the FG slave node ",nodename
2090         endif
2091 #endif
2092         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2093         lenpre = lentmp+lenpre+1
2094       endif
2095       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2096 C Get the names and open the input files
2097 #if defined(WINIFL) || defined(WINPGI)
2098       open(1,file=pref_orig(:ilen(pref_orig))//
2099      &  '.inp',status='old',readonly,shared)
2100        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2101 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2102 C Get parameter filenames and open the parameter files.
2103       call getenv_loc('BONDPAR',bondname)
2104       open (ibond,file=bondname,status='old',readonly,shared)
2105       call getenv_loc('THETPAR',thetname)
2106       open (ithep,file=thetname,status='old',readonly,shared)
2107       call getenv_loc('ROTPAR',rotname)
2108       open (irotam,file=rotname,status='old',readonly,shared)
2109       call getenv_loc('TORPAR',torname)
2110       open (itorp,file=torname,status='old',readonly,shared)
2111       call getenv_loc('TORDPAR',tordname)
2112       open (itordp,file=tordname,status='old',readonly,shared)
2113       call getenv_loc('FOURIER',fouriername)
2114       open (ifourier,file=fouriername,status='old',readonly,shared)
2115       call getenv_loc('ELEPAR',elename)
2116       open (ielep,file=elename,status='old',readonly,shared)
2117       call getenv_loc('SIDEPAR',sidename)
2118       open (isidep,file=sidename,status='old',readonly,shared)
2119       call getenv_loc('LIPTRANPAR',liptranname)
2120       open (iliptranpar,file=liptranname,status='old',readonly,shared)
2121 #elif (defined CRAY) || (defined AIX)
2122       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2123      &  action='read')
2124 c      print *,"Processor",myrank," opened file 1"
2125       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2126 c      print *,"Processor",myrank," opened file 9"
2127 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2128 C Get parameter filenames and open the parameter files.
2129       call getenv_loc('BONDPAR',bondname)
2130       open (ibond,file=bondname,status='old',action='read')
2131 c      print *,"Processor",myrank," opened file IBOND"
2132       call getenv_loc('THETPAR',thetname)
2133       open (ithep,file=thetname,status='old',action='read')
2134 #ifndef CRYST_THETA
2135       call getenv_loc('THETPARPDB',thetname_pdb)
2136       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2137 #endif
2138 c      print *,"Processor",myrank," opened file ITHEP"
2139       call getenv_loc('ROTPAR',rotname)
2140       open (irotam,file=rotname,status='old',action='read')
2141 #ifndef CRYST_SC
2142       call getenv_loc('ROTPARPDB',rotname_pdb)
2143       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2144 #endif
2145 c      print *,"Processor",myrank," opened file IROTAM"
2146       call getenv_loc('TORPAR',torname)
2147       open (itorp,file=torname,status='old',action='read')
2148 c      print *,"Processor",myrank," opened file ITORP"
2149       call getenv_loc('TORDPAR',tordname)
2150       open (itordp,file=tordname,status='old',action='read')
2151 c      print *,"Processor",myrank," opened file ITORDP"
2152       call getenv_loc('SCCORPAR',sccorname)
2153       open (isccor,file=sccorname,status='old',action='read')
2154 c      print *,"Processor",myrank," opened file ISCCOR"
2155       call getenv_loc('FOURIER',fouriername)
2156       open (ifourier,file=fouriername,status='old',action='read')
2157 c      print *,"Processor",myrank," opened file IFOURIER"
2158       call getenv_loc('ELEPAR',elename)
2159       open (ielep,file=elename,status='old',action='read')
2160 c      print *,"Processor",myrank," opened file IELEP"
2161       call getenv_loc('SIDEPAR',sidename)
2162       open (isidep,file=sidename,status='old',action='read')
2163       call getenv_loc('LIPTRANPAR',liptranname)
2164       open (iliptranpar,file=liptranname,status='old',action='read')
2165 c      print *,"Processor",myrank," opened file ISIDEP"
2166 c      print *,"Processor",myrank," opened parameter files"
2167 #elif (defined G77)
2168       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2169       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2170 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2171 C Get parameter filenames and open the parameter files.
2172       call getenv_loc('BONDPAR',bondname)
2173       open (ibond,file=bondname,status='old')
2174       call getenv_loc('THETPAR',thetname)
2175       open (ithep,file=thetname,status='old')
2176 #ifndef CRYST_THETA
2177       call getenv_loc('THETPARPDB',thetname_pdb)
2178       open (ithep_pdb,file=thetname_pdb,status='old')
2179 #endif
2180       call getenv_loc('ROTPAR',rotname)
2181       open (irotam,file=rotname,status='old')
2182 #ifndef CRYST_SC
2183       call getenv_loc('ROTPARPDB',rotname_pdb)
2184       open (irotam_pdb,file=rotname_pdb,status='old')
2185 #endif
2186       call getenv_loc('TORPAR',torname)
2187       open (itorp,file=torname,status='old')
2188       call getenv_loc('TORDPAR',tordname)
2189       open (itordp,file=tordname,status='old')
2190       call getenv_loc('SCCORPAR',sccorname)
2191       open (isccor,file=sccorname,status='old')
2192       call getenv_loc('FOURIER',fouriername)
2193       open (ifourier,file=fouriername,status='old')
2194       call getenv_loc('ELEPAR',elename)
2195       open (ielep,file=elename,status='old')
2196       call getenv_loc('SIDEPAR',sidename)
2197       open (isidep,file=sidename,status='old')
2198       call getenv_loc('LIPTRANPAR',liptranname)
2199       open (iliptranpar,file=liptranname,status='old')
2200 #else
2201       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2202      &  readonly)
2203        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2204 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2205 C Get parameter filenames and open the parameter files.
2206       call getenv_loc('BONDPAR',bondname)
2207       open (ibond,file=bondname,status='old',readonly)
2208       call getenv_loc('THETPAR',thetname)
2209       open (ithep,file=thetname,status='old',readonly)
2210       call getenv_loc('ROTPAR',rotname)
2211       open (irotam,file=rotname,status='old',readonly)
2212       call getenv_loc('TORPAR',torname)
2213       open (itorp,file=torname,status='old',readonly)
2214       call getenv_loc('TORDPAR',tordname)
2215       open (itordp,file=tordname,status='old',readonly)
2216       call getenv_loc('SCCORPAR',sccorname)
2217       open (isccor,file=sccorname,status='old',readonly)
2218 #ifndef CRYST_THETA
2219       call getenv_loc('THETPARPDB',thetname_pdb)
2220       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2221 #endif
2222       call getenv_loc('FOURIER',fouriername)
2223       open (ifourier,file=fouriername,status='old',readonly)
2224       call getenv_loc('ELEPAR',elename)
2225       open (ielep,file=elename,status='old',readonly)
2226       call getenv_loc('SIDEPAR',sidename)
2227       open (isidep,file=sidename,status='old',readonly)
2228       call getenv_loc('LIPTRANPAR',liptranname)
2229       open (iliptranpar,file=liptranname,status='old',action='read')
2230 #ifndef CRYST_SC
2231       call getenv_loc('ROTPARPDB',rotname_pdb)
2232       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2233 #endif
2234 #endif
2235 #ifdef TUBE
2236       call getenv_loc('TUBEPAR',tubename)
2237 #if defined(WINIFL) || defined(WINPGI)
2238       open (itube,file=tubename,status='old',readonly,shared)
2239 #elif (defined CRAY)  || (defined AIX)
2240       open (itube,file=tubename,status='old',action='read')
2241 #elif (defined G77)
2242       open (itube,file=tubename,status='old')
2243 #else
2244       open (itube,file=tubename,status='old',readonly)
2245 #endif
2246 #endif
2247 #ifndef OLDSCP
2248 C
2249 C 8/9/01 In the newest version SCp interaction constants are read from a file
2250 C Use -DOLDSCP to use hard-coded constants instead.
2251 C
2252       call getenv_loc('SCPPAR',scpname)
2253 #if defined(WINIFL) || defined(WINPGI)
2254       open (iscpp,file=scpname,status='old',readonly,shared)
2255 #elif (defined CRAY)  || (defined AIX)
2256       open (iscpp,file=scpname,status='old',action='read')
2257 #elif (defined G77)
2258       open (iscpp,file=scpname,status='old')
2259 #else
2260       open (iscpp,file=scpname,status='old',readonly)
2261 #endif
2262 #endif
2263       call getenv_loc('PATTERN',patname)
2264 #if defined(WINIFL) || defined(WINPGI)
2265       open (icbase,file=patname,status='old',readonly,shared)
2266 #elif (defined CRAY)  || (defined AIX)
2267       open (icbase,file=patname,status='old',action='read')
2268 #elif (defined G77)
2269       open (icbase,file=patname,status='old')
2270 #else
2271       open (icbase,file=patname,status='old',readonly)
2272 #endif
2273 #ifdef MPI
2274 C Open output file only for CG processes
2275 c      print *,"Processor",myrank," fg_rank",fg_rank
2276       if (fg_rank.eq.0) then
2277
2278       if (nodes.eq.1) then
2279         npos=3
2280       else
2281         npos = dlog10(dfloat(nodes-1))+1
2282       endif
2283       if (npos.lt.3) npos=3
2284       write (liczba,'(i1)') npos
2285       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2286      &  //')'
2287       write (liczba,form) me
2288       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2289      &  liczba(:ilen(liczba))
2290       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2291      &  //'.int'
2292       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2293      &  //'.pdb'
2294       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2295      &  liczba(:ilen(liczba))//'.mol2'
2296       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2297      &  liczba(:ilen(liczba))//'.stat'
2298       if (lentmp.gt.0)
2299      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2300      &      //liczba(:ilen(liczba))//'.stat')
2301       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2302      &  //'.rst'
2303       if(usampl) then
2304           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2305      & liczba(:ilen(liczba))//'.const'
2306       endif
2307
2308       endif
2309 #else
2310       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2311       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2312       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2313       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2314       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2315       if (lentmp.gt.0)
2316      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2317      &    //'.stat')
2318       rest2name=prefix(:ilen(prefix))//'.rst'
2319       if(usampl) then
2320          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2321       endif
2322 #endif
2323 #if defined(AIX) || defined(PGI) || defined(CRAY)
2324       if (me.eq.king .or. .not. out1file) then
2325          open(iout,file=outname,status='unknown')
2326       else
2327          open(iout,file="/dev/null",status="unknown")
2328       endif
2329 c#define DEBUG
2330 #ifdef DEBUG
2331       if (fg_rank.gt.0) then
2332         write (liczba,'(i3.3)') myrank/nfgtasks
2333         write (ll,'(bz,i3.3)') fg_rank
2334         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2335      &   status='unknown')
2336       endif
2337 #endif
2338 c#undef DEBUG
2339       if(me.eq.king) then
2340        open(igeom,file=intname,status='unknown',position='append')
2341        open(ipdb,file=pdbname,status='unknown')
2342        open(imol2,file=mol2name,status='unknown')
2343        open(istat,file=statname,status='unknown',position='append')
2344       else
2345 c1out       open(iout,file=outname,status='unknown')
2346       endif
2347 #else
2348       if (me.eq.king .or. .not.out1file)
2349      &    open(iout,file=outname,status='unknown')
2350 #ifdef DEBUG
2351       if (fg_rank.gt.0) then
2352         write (liczba,'(i3.3)') myrank/nfgtasks
2353         write (ll,'(bz,i3.3)') fg_rank
2354         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2355      &   status='unknown')
2356       endif
2357 #endif
2358       if(me.eq.king) then
2359        open(igeom,file=intname,status='unknown',access='append')
2360        open(ipdb,file=pdbname,status='unknown')
2361        open(imol2,file=mol2name,status='unknown')
2362        open(istat,file=statname,status='unknown',access='append')
2363       else
2364 c1out       open(iout,file=outname,status='unknown')
2365       endif
2366 #endif
2367       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2368       csa_seed=prefix(:lenpre)//'.CSA.seed'
2369       csa_history=prefix(:lenpre)//'.CSA.history'
2370       csa_bank=prefix(:lenpre)//'.CSA.bank'
2371       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2372       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2373       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2374 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2375       csa_int=prefix(:lenpre)//'.int'
2376       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2377       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2378       csa_in=prefix(:lenpre)//'.CSA.in'
2379 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2380 C Write file names
2381       if (me.eq.king)then
2382       write (iout,'(80(1h-))')
2383       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2384       write (iout,'(80(1h-))')
2385       write (iout,*) "Input file                      : ",
2386      &  pref_orig(:ilen(pref_orig))//'.inp'
2387       write (iout,*) "Output file                     : ",
2388      &  outname(:ilen(outname))
2389       write (iout,*)
2390       write (iout,*) "Sidechain potential file        : ",
2391      &  sidename(:ilen(sidename))
2392 #ifndef OLDSCP
2393       write (iout,*) "SCp potential file              : ",
2394      &  scpname(:ilen(scpname))
2395 #endif
2396       write (iout,*) "Electrostatic potential file    : ",
2397      &  elename(:ilen(elename))
2398       write (iout,*) "Cumulant coefficient file       : ",
2399      &  fouriername(:ilen(fouriername))
2400       write (iout,*) "Torsional parameter file        : ",
2401      &  torname(:ilen(torname))
2402       write (iout,*) "Double torsional parameter file : ",
2403      &  tordname(:ilen(tordname))
2404       write (iout,*) "SCCOR parameter file : ",
2405      &  sccorname(:ilen(sccorname))
2406       write (iout,*) "Bond & inertia constant file    : ",
2407      &  bondname(:ilen(bondname))
2408       write (iout,*) "Bending-torsion parameter file  : ",
2409      &  thetname(:ilen(thetname))
2410       write (iout,*) "Rotamer parameter file          : ",
2411      &  rotname(:ilen(rotname))
2412       write (iout,*) "Thetpdb parameter file          : ",
2413      &  thetname_pdb(:ilen(thetname_pdb))
2414       write (iout,*) "Threading database              : ",
2415      &  patname(:ilen(patname))
2416       if (lentmp.ne.0)
2417      &write (iout,*)" DIRTMP                          : ",
2418      &  tmpdir(:lentmp)
2419       write (iout,'(80(1h-))')
2420       endif
2421       return
2422       end
2423 c----------------------------------------------------------------------------
2424       subroutine card_concat(card)
2425       implicit real*8 (a-h,o-z)
2426       include 'DIMENSIONS'
2427       include 'COMMON.IOUNITS'
2428       character*(*) card
2429       character*80 karta,ucase
2430       external ilen
2431       read (inp,'(a)') karta
2432       karta=ucase(karta)
2433       card=' '
2434       do while (karta(80:80).eq.'&')
2435         card=card(:ilen(card)+1)//karta(:79)
2436         read (inp,'(a)') karta
2437         karta=ucase(karta)
2438       enddo
2439       card=card(:ilen(card)+1)//karta
2440       return
2441       end
2442 c----------------------------------------------------------------------------------
2443       subroutine readrst
2444       implicit real*8 (a-h,o-z)
2445       include 'DIMENSIONS'
2446       include 'COMMON.CHAIN'
2447       include 'COMMON.IOUNITS'
2448       include 'COMMON.MD'
2449       open(irest2,file=rest2name,status='unknown')
2450       read(irest2,*) totT,EK,potE,totE,t_bath
2451       totTafm=totT
2452       do i=1,2*nres
2453          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2454       enddo
2455       do i=1,2*nres
2456          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2457       enddo
2458       if(usampl) then
2459              read (irest2,*) iset
2460       endif
2461       close(irest2)
2462       return
2463       end
2464 c---------------------------------------------------------------------------------
2465       subroutine read_fragments
2466       implicit real*8 (a-h,o-z)
2467       include 'DIMENSIONS'
2468 #ifdef MPI
2469       include 'mpif.h'
2470 #endif
2471       include 'COMMON.SETUP'
2472       include 'COMMON.CHAIN'
2473       include 'COMMON.IOUNITS'
2474       include 'COMMON.MD'
2475       include 'COMMON.CONTROL'
2476       read(inp,*) nset,nfrag,npair,nfrag_back
2477       loc_qlike=(nfrag_back.lt.0)
2478       nfrag_back=iabs(nfrag_back)
2479       if(me.eq.king.or..not.out1file)
2480      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2481      &  " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2482       do iset=1,nset
2483          read(inp,*) mset(iset)
2484        do i=1,nfrag
2485          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2486      &     qinfrag(i,iset)
2487          if(me.eq.king.or..not.out1file)
2488      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2489      &     ifrag(2,i,iset), qinfrag(i,iset)
2490        enddo
2491        do i=1,npair
2492         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2493      &    qinpair(i,iset)
2494         if(me.eq.king.or..not.out1file)
2495      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2496      &    ipair(2,i,iset), qinpair(i,iset)
2497        enddo
2498        if (loc_qlike) then
2499        do i=1,nfrag_back
2500         read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2501      &     wfrag_back(2,i,iset),qin_back(2,i,iset),
2502      &     wfrag_back(3,i,iset),qin_back(3,i,iset),
2503      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2504         if(me.eq.king.or..not.out1file)
2505      &   write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2506      &   wfrag_back(2,i,iset),qin_back(3,i,iset),
2507      &   wfrag_back(3,i,iset),qin_back(3,i,iset),
2508      &   ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2509        enddo
2510        else
2511        do i=1,nfrag_back
2512         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2513      &     wfrag_back(3,i,iset),
2514      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2515         if(me.eq.king.or..not.out1file)
2516      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2517      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2518        enddo
2519        endif
2520       enddo
2521       return
2522       end
2523 C---------------------------------------------------------------------------
2524       subroutine read_afminp
2525             implicit real*8 (a-h,o-z)
2526       include 'DIMENSIONS'
2527 #ifdef MPI
2528       include 'mpif.h'
2529 #endif
2530       include 'COMMON.SETUP'
2531       include 'COMMON.CONTROL'
2532       include 'COMMON.CHAIN'
2533       include 'COMMON.IOUNITS'
2534       include 'COMMON.SBRIDGE'
2535       character*320 afmcard
2536       print *, "wchodze"
2537       call card_concat(afmcard)
2538       call readi(afmcard,"BEG",afmbeg,0)
2539       call readi(afmcard,"END",afmend,0)
2540       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2541       call reada(afmcard,"VEL",velAFMconst,0.0d0)
2542       print *,'FORCE=' ,forceAFMconst
2543 CCCC NOW PROPERTIES FOR AFM
2544        distafminit=0.0d0
2545        do i=1,3
2546         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2547        enddo
2548         distafminit=dsqrt(distafminit)
2549         print *,'initdist',distafminit
2550       return
2551       end
2552 c-------------------------------------------------------------------------------
2553       subroutine read_saxs_constr
2554       implicit real*8 (a-h,o-z)
2555       include 'DIMENSIONS'
2556 #ifdef MPI
2557       include 'mpif.h'
2558 #endif
2559       include 'COMMON.SETUP'
2560       include 'COMMON.CONTROL'
2561       include 'COMMON.CHAIN'
2562       include 'COMMON.IOUNITS'
2563       include 'COMMON.SBRIDGE'
2564       double precision cm(3)
2565 c      read(inp,*) nsaxs
2566       write (iout,*) "Calling read_saxs nsaxs",nsaxs
2567       call flush(iout)
2568       if (saxs_mode.eq.0) then
2569 c SAXS distance distribution
2570       do i=1,nsaxs
2571         read(inp,*) distsaxs(i),Psaxs(i)
2572       enddo
2573       Cnorm = 0.0d0
2574       do i=1,nsaxs
2575         Cnorm = Cnorm + Psaxs(i)
2576       enddo
2577       write (iout,*) "Cnorm",Cnorm
2578       do i=1,nsaxs
2579         Psaxs(i)=Psaxs(i)/Cnorm
2580       enddo
2581       write (iout,*) "Normalized distance distribution from SAXS"
2582       do i=1,nsaxs
2583         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2584       enddo
2585       Wsaxs0=0.0d0
2586       do i=1,nsaxs
2587         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2588       enddo
2589       write (iout,*) "Wsaxs0",Wsaxs0
2590       else
2591 c SAXS "spheres".
2592       do i=1,nsaxs
2593         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2594       enddo
2595       do j=1,3
2596         cm(j)=0.0d0
2597       enddo
2598       do i=1,nsaxs
2599         do j=1,3
2600           cm(j)=cm(j)+Csaxs(j,i)
2601         enddo
2602       enddo
2603       do j=1,3
2604         cm(j)=cm(j)/nsaxs
2605       enddo
2606       do i=1,nsaxs
2607         do j=1,3
2608           Csaxs(j,i)=Csaxs(j,i)-cm(j)
2609         enddo
2610       enddo
2611       write (iout,*) "SAXS sphere coordinates"
2612       do i=1,nsaxs
2613         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2614       enddo
2615       endif
2616       return
2617       end
2618
2619 c-------------------------------------------------------------------------------
2620       subroutine read_dist_constr
2621       implicit real*8 (a-h,o-z)
2622       include 'DIMENSIONS'
2623 #ifdef MPI
2624       include 'mpif.h'
2625 #endif
2626       include 'COMMON.SETUP'
2627       include 'COMMON.CONTROL'
2628       include 'COMMON.CHAIN'
2629       include 'COMMON.IOUNITS'
2630       include 'COMMON.SBRIDGE'
2631       include 'COMMON.INTERACT'
2632       integer ifrag_(2,100),ipair_(2,1000)
2633       double precision wfrag_(100),wpair_(1000)
2634       character*5000 controlcard
2635       logical normalize,next
2636       integer restr_type
2637       double precision scal_bfac
2638       double precision xlink(4,0:4) /
2639 c           a          b       c     sigma
2640      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
2641      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
2642      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
2643      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
2644      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
2645 c      print *, "WCHODZE"
2646       write (iout,*) "Calling read_dist_constr"
2647 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2648 c      call flush(iout)
2649       restr_on_coord=.false.
2650       next=.true.
2651
2652       npeak=0
2653       ipeak=0
2654       nhpb_peak=0
2655
2656       DO WHILE (next)
2657
2658       call card_concat(controlcard)
2659       next = index(controlcard,"NEXT").gt.0
2660       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2661       write (iout,*) "restr_type",restr_type
2662       call readi(controlcard,"NFRAG",nfrag_,0)
2663       call readi(controlcard,"NFRAG",nfrag_,0)
2664       call readi(controlcard,"NPAIR",npair_,0)
2665       call readi(controlcard,"NDIST",ndist_,0)
2666       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2667       call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2668       if (restr_type.eq.10)
2669      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2670       if (restr_type.eq.12)
2671      &  call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2672       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2673       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2674       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2675       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2676       normalize = index(controlcard,"NORMALIZE").gt.0
2677       write (iout,*) "WBOLTZD",wboltzd
2678       write (iout,*) "SCAL_PEAK",scal_peak
2679       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2680       write (iout,*) "IFRAG"
2681       do i=1,nfrag_
2682         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2683       enddo
2684       write (iout,*) "IPAIR"
2685       do i=1,npair_
2686         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2687       enddo
2688       if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2689      & write (iout,*)
2690      &   "Distance restraints as generated from reference structure"
2691       do i=1,nfrag_
2692         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2693         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2694      &    ifrag_(2,i)=nstart_sup+nsup-1
2695 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2696 c        call flush(iout)
2697         if (wfrag_(i).eq.0.0d0) cycle
2698         do j=ifrag_(1,i),ifrag_(2,i)-1
2699           do k=j+1,ifrag_(2,i)
2700 c            write (iout,*) "j",j," k",k
2701             ddjk=dist(j,k)
2702             if (restr_type.eq.1) then
2703               nhpb=nhpb+1
2704               irestr_type(nhpb)=1
2705               ihpb(nhpb)=j
2706               jhpb(nhpb)=k
2707               dhpb(nhpb)=ddjk
2708               forcon(nhpb)=wfrag_(i)
2709             else if (constr_dist.eq.2) then
2710               if (ddjk.le.dist_cut) then
2711                 nhpb=nhpb+1
2712                 irestr_type(nhpb)=1
2713                 ihpb(nhpb)=j
2714                 jhpb(nhpb)=k
2715                 dhpb(nhpb)=ddjk
2716                 forcon(nhpb)=wfrag_(i)
2717               endif
2718             else if (restr_type.eq.3) then
2719               nhpb=nhpb+1
2720               irestr_type(nhpb)=1
2721               ihpb(nhpb)=j
2722               jhpb(nhpb)=k
2723               dhpb(nhpb)=ddjk
2724               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2725             endif
2726 #ifdef MPI
2727             if (.not.out1file .or. me.eq.king)
2728      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2729      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2730 #else
2731             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2732      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2733 #endif
2734           enddo
2735         enddo
2736       enddo
2737       do i=1,npair_
2738         if (wpair_(i).eq.0.0d0) cycle
2739         ii = ipair_(1,i)
2740         jj = ipair_(2,i)
2741         if (ii.gt.jj) then
2742           itemp=ii
2743           ii=jj
2744           jj=itemp
2745         endif
2746         do j=ifrag_(1,ii),ifrag_(2,ii)
2747           do k=ifrag_(1,jj),ifrag_(2,jj)
2748             ddjk=dist(j,k)
2749             if (restr_type.eq.1) then
2750               nhpb=nhpb+1
2751               irestr_type(nhpb)=1
2752               ihpb(nhpb)=j
2753               jhpb(nhpb)=k
2754               dhpb(nhpb)=ddjk
2755               forcon(nhpb)=wpair_(i)
2756             else if (constr_dist.eq.2) then
2757               if (ddjk.le.dist_cut) then
2758                 nhpb=nhpb+1
2759                 irestr_type(nhpb)=1
2760                 ihpb(nhpb)=j
2761                 jhpb(nhpb)=k
2762                 dhpb(nhpb)=ddjk
2763                 forcon(nhpb)=wpair_(i)
2764               endif
2765             else if (restr_type.eq.3) then
2766               nhpb=nhpb+1
2767               irestr_type(nhpb)=1
2768               ihpb(nhpb)=j
2769               jhpb(nhpb)=k
2770               dhpb(nhpb)=ddjk
2771               forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2772             endif
2773 #ifdef MPI
2774             if (.not.out1file .or. me.eq.king)
2775      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2776      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2777 #else
2778             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2779      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2780 #endif
2781           enddo
2782         enddo
2783       enddo
2784
2785 c      print *,ndist_
2786       write (iout,*) "Distance restraints as read from input"
2787       do i=1,ndist_
2788         if (restr_type.eq.12) then
2789           read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2790      &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2791      &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2792      &    fordepth_peak(nhpb_peak+1),npeak
2793 c          write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2794 c     &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2795 c     &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2796 c     &    fordepth_peak(nhpb_peak+1),npeak
2797           if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2798      &      fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2799           nhpb_peak=nhpb_peak+1
2800           irestr_type_peak(nhpb_peak)=12
2801           if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2802           ipeak(2,npeak)=i
2803 #ifdef MPI
2804           if (.not.out1file .or. me.eq.king)
2805      &    write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2806      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2807      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2808      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2809      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2810 #else
2811           write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2812      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2813      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2814      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2815      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2816 #endif
2817           if (ibecarb_peak(nhpb_peak).eq.3) then
2818             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2819           else if (ibecarb_peak(nhpb_peak).eq.2) then
2820             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2821           else if (ibecarb_peak(nhpb_peak).eq.1) then
2822             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2823             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2824           endif
2825         else if (restr_type.eq.11) then
2826           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2827      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2828 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2829           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2830           nhpb=nhpb+1
2831           irestr_type(nhpb)=11
2832 #ifdef MPI
2833           if (.not.out1file .or. me.eq.king)
2834      &    write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2835      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2836      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2837 #else
2838           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2839      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2840      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2841 #endif
2842 c          if (ibecarb(nhpb).gt.0) then
2843 c            ihpb(nhpb)=ihpb(nhpb)+nres
2844 c            jhpb(nhpb)=jhpb(nhpb)+nres
2845 c          endif
2846          if (ibecarb(nhpb).eq.3) then
2847             ihpb(nhpb)=ihpb(nhpb)+nres
2848           else if (ibecarb(nhpb).eq.2) then
2849             ihpb(nhpb)=ihpb(nhpb)+nres
2850           else if (ibecarb(nhpb).eq.1) then
2851             ihpb(nhpb)=ihpb(nhpb)+nres
2852             jhpb(nhpb)=jhpb(nhpb)+nres
2853           endif
2854         else if (restr_type.eq.10) then
2855 c Cross-lonk Markov-like potential
2856           call card_concat(controlcard)
2857           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2858           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2859           ibecarb(nhpb+1)=0
2860           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2861           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2862           if (index(controlcard,"ZL").gt.0) then
2863             link_type=1
2864           else if (index(controlcard,"ADH").gt.0) then
2865             link_type=2
2866           else if (index(controlcard,"PDH").gt.0) then
2867             link_type=3
2868           else if (index(controlcard,"DSS").gt.0) then
2869             link_type=4
2870           else
2871             link_type=0
2872           endif
2873           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2874      &       xlink(1,link_type))
2875           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2876      &       xlink(2,link_type))
2877           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2878      &       xlink(3,link_type))
2879           call reada(controlcard,"SIGMA",forcon(nhpb+1),
2880      &       xlink(4,link_type))
2881           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2882 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2883 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2884           if (forcon(nhpb+1).le.0.0d0 .or.
2885      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2886           nhpb=nhpb+1
2887           irestr_type(nhpb)=10
2888           if (ibecarb(nhpb).eq.3) then
2889             jhpb(nhpb)=jhpb(nhpb)+nres
2890           else if (ibecarb(nhpb).eq.2) then
2891             ihpb(nhpb)=ihpb(nhpb)+nres
2892           else if (ibecarb(nhpb).eq.1) then
2893             ihpb(nhpb)=ihpb(nhpb)+nres
2894             jhpb(nhpb)=jhpb(nhpb)+nres
2895           endif
2896 #ifdef MPI
2897           if (.not.out1file .or. me.eq.king)
2898      &    write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2899      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2900      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2901      &     irestr_type(nhpb)
2902 #else
2903           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2904      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2905      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2906      &     irestr_type(nhpb)
2907 #endif
2908         else
2909 C        print *,"in else"
2910           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2911      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2912           if (forcon(nhpb+1).gt.0.0d0) then
2913           nhpb=nhpb+1
2914           if (dhpb1(nhpb).eq.0.0d0) then
2915             irestr_type(nhpb)=1
2916           else
2917             irestr_type(nhpb)=2
2918           endif
2919           if (ibecarb(nhpb).eq.3) then
2920             jhpb(nhpb)=jhpb(nhpb)+nres
2921           else if (ibecarb(nhpb).eq.2) then
2922             ihpb(nhpb)=ihpb(nhpb)+nres
2923           else if (ibecarb(nhpb).eq.1) then
2924             ihpb(nhpb)=ihpb(nhpb)+nres
2925             jhpb(nhpb)=jhpb(nhpb)+nres
2926           endif
2927           if (dhpb(nhpb).eq.0.0d0)
2928      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2929         endif
2930 #ifdef MPI
2931           if (.not.out1file .or. me.eq.king)
2932      &    write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2933      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2934 #else
2935           write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2936      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2937 #endif
2938         endif
2939 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2940 C        if (forcon(nhpb+1).gt.0.0d0) then
2941 C          nhpb=nhpb+1
2942 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2943       enddo
2944
2945       if (restr_type.eq.4) then
2946         write (iout,*) "The BFAC array"
2947         do i=nnt,nct
2948           write (iout,'(i5,f10.5)') i,bfac(i)
2949         enddo
2950         do i=nnt,nct
2951           if (itype(i).eq.ntyp1) cycle
2952           do j=nnt,i-1
2953             if (itype(j).eq.ntyp1) cycle
2954             if (itype(i).eq.10) then
2955               iiend=0
2956             else
2957               iiend=1
2958             endif
2959             if (itype(j).eq.10) then
2960               jjend=0
2961             else
2962               jjend=1
2963             endif
2964             kk=0
2965             do ii=0,iiend
2966             do jj=0,jjend
2967             nhpb=nhpb+1
2968             irestr_type(nhpb)=1
2969             forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2970             irestr_type(nhpb)=1
2971             ibecarb(nhpb)=kk
2972             if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2973             ihpb(nhpb)=i+nres*ii
2974             jhpb(nhpb)=j+nres*jj
2975             dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2976 #ifdef MPI
2977             if (.not.out1file .or. me.eq.king) then
2978             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2979      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2980      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2981      &       irestr_type(nhpb)
2982             endif
2983 #else
2984             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2985      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2986      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2987      &       irestr_type(nhpb)
2988 #endif
2989             kk=kk+1
2990           enddo
2991           enddo
2992           enddo
2993         enddo
2994       endif
2995
2996       if (restr_type.eq.5) then
2997         restr_on_coord=.true.
2998         do i=nnt,nct
2999           if (itype(i).eq.ntyp1) cycle
3000           bfac(i)=(scal_bfac/bfac(i))**2
3001         enddo
3002       endif
3003
3004       ENDDO ! next
3005
3006       fordepthmax=0.0d0
3007       if (normalize) then
3008         do i=nss+1,nhpb
3009           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
3010      &      fordepthmax=fordepth(i)
3011         enddo
3012         do i=nss+1,nhpb
3013           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
3014         enddo
3015       endif
3016       return
3017       end
3018 c-------------------------------------------------------------------------------
3019       subroutine read_constr_homology
3020
3021       include 'DIMENSIONS'
3022 #ifdef MPI
3023       include 'mpif.h'
3024 #endif
3025       include 'COMMON.SETUP'
3026       include 'COMMON.CONTROL'
3027       include 'COMMON.CHAIN'
3028       include 'COMMON.IOUNITS'
3029       include 'COMMON.MD'
3030       include 'COMMON.GEO'
3031       include 'COMMON.INTERACT'
3032       include 'COMMON.NAMES'
3033 c
3034 c For new homol impl
3035 c
3036       include 'COMMON.VAR'
3037 c
3038
3039 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
3040 c    &                 dist_cut
3041 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
3042 c    &    sigma_odl_temp(maxres,maxres,max_template)
3043       character*2 kic2
3044       character*24 model_ki_dist, model_ki_angle
3045       character*500 controlcard
3046       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
3047       integer ilen
3048       external ilen
3049       logical liiflag
3050 c
3051 c     FP - Nov. 2014 Temporary specifications for new vars
3052 c
3053       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
3054      &                 rescore3_tmp
3055       double precision, dimension (max_template,maxres) :: rescore
3056       double precision, dimension (max_template,maxres) :: rescore2
3057       double precision, dimension (max_template,maxres) :: rescore3
3058       character*24 pdbfile,tpl_k_rescore
3059 c -----------------------------------------------------------------
3060 c Reading multiple PDB ref structures and calculation of retraints
3061 c not using pre-computed ones stored in files model_ki_{dist,angle}
3062 c FP (Nov., 2014)
3063 c -----------------------------------------------------------------
3064 c
3065 c
3066 c Alternative: reading from input
3067       call card_concat(controlcard)
3068       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
3069       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
3070       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
3071       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
3072       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3073       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3074       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3075       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3076       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3077       if(.not.read2sigma.and.start_from_model) then
3078           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3079      &      write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3080           start_from_model=.false.
3081       endif
3082       if(start_from_model .and. (me.eq.king .or. .not. out1file))
3083      &    write(iout,*) 'START_FROM_MODELS is ON'
3084       if(start_from_model .and. rest) then
3085         if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3086          write(iout,*) 'START_FROM_MODELS is OFF'
3087          write(iout,*) 'remove restart keyword from input'
3088         endif
3089       endif
3090       if (homol_nset.gt.1)then
3091          call card_concat(controlcard)
3092          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3093          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3094           write(iout,*) "iset homology_weight "
3095           do i=1,homol_nset
3096            write(iout,*) i,waga_homology(i)
3097           enddo
3098          endif
3099          iset=mod(kolor,homol_nset)+1
3100       else
3101        iset=1
3102        waga_homology(1)=1.0
3103       endif
3104
3105 cd      write (iout,*) "nnt",nnt," nct",nct
3106 cd      call flush(iout)
3107
3108
3109       lim_odl=0
3110       lim_dih=0
3111 c
3112 c      write(iout,*) 'nnt=',nnt,'nct=',nct
3113 c
3114       do i = nnt,nct
3115         do k=1,constr_homology
3116          idomain(k,i)=0
3117         enddo
3118       enddo
3119
3120       ii=0
3121       do i = nnt,nct-2
3122         do j=i+2,nct
3123         ii=ii+1
3124         ii_in_use(ii)=0
3125         enddo
3126       enddo
3127
3128       if (read_homol_frag) then
3129         call read_klapaucjusz
3130       else
3131
3132       do k=1,constr_homology
3133
3134         read(inp,'(a)') pdbfile
3135         if(me.eq.king .or. .not. out1file)
3136      &    write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3137      &  pdbfile(:ilen(pdbfile))
3138         open(ipdbin,file=pdbfile,status='old',err=33)
3139         goto 34
3140   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3141      &  pdbfile(:ilen(pdbfile))
3142         stop
3143   34    continue
3144 c        print *,'Begin reading pdb data'
3145 c
3146 c Files containing res sim or local scores (former containing sigmas)
3147 c
3148
3149         write(kic2,'(bz,i2.2)') k
3150
3151         tpl_k_rescore="template"//kic2//".sco"
3152
3153         unres_pdb=.false.
3154         if (read2sigma) then
3155           call readpdb_template(k)
3156         else
3157           call readpdb
3158         endif
3159 c
3160 c     Distance restraints
3161 c
3162 c          ... --> odl(k,ii)
3163 C Copy the coordinates from reference coordinates (?)
3164         do i=1,2*nres
3165           do j=1,3
3166             c(j,i)=cref(j,i)
3167 c           write (iout,*) "c(",j,i,") =",c(j,i)
3168           enddo
3169         enddo
3170 c
3171 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3172 c
3173 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3174           open (ientin,file=tpl_k_rescore,status='old')
3175           if (nnt.gt.1) rescore(k,1)=0.0d0
3176           do irec=nnt,nct ! loop for reading res sim
3177             if (read2sigma) then
3178              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3179      &                                rescore3_tmp,idomain_tmp
3180              i_tmp=i_tmp+nnt-1
3181              idomain(k,i_tmp)=idomain_tmp
3182              rescore(k,i_tmp)=rescore_tmp
3183              rescore2(k,i_tmp)=rescore2_tmp
3184              rescore3(k,i_tmp)=rescore3_tmp
3185              if (.not. out1file .or. me.eq.king)
3186      &        write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3187      &                      i_tmp,rescore2_tmp,rescore_tmp,
3188      &                                rescore3_tmp,idomain_tmp
3189             else
3190              idomain(k,irec)=1
3191              read (ientin,*,end=1401) rescore_tmp
3192
3193 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3194              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3195 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3196             endif
3197           enddo
3198  1401   continue
3199         close (ientin)
3200         if (waga_dist.ne.0.0d0) then
3201           ii=0
3202           do i = nnt,nct-2
3203             do j=i+2,nct
3204
3205               x12=c(1,i)-c(1,j)
3206               y12=c(2,i)-c(2,j)
3207               z12=c(3,i)-c(3,j)
3208               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3209 c              write (iout,*) k,i,j,distal,dist2_cut
3210
3211             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3212      &            .and. distal.le.dist2_cut ) then
3213
3214               ii=ii+1
3215               ii_in_use(ii)=1
3216               l_homo(k,ii)=.true.
3217
3218 c             write (iout,*) "k",k
3219 c             write (iout,*) "i",i," j",j," constr_homology",
3220 c    &                       constr_homology
3221               ires_homo(ii)=i
3222               jres_homo(ii)=j
3223               odl(k,ii)=distal
3224               if (read2sigma) then
3225                 sigma_odl(k,ii)=0
3226                 do ik=i,j
3227                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3228                 enddo
3229                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3230                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3231      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3232               else
3233                 if (odl(k,ii).le.dist_cut) then
3234                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3235                 else
3236 #ifdef OLDSIGMA
3237                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3238      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3239 #else
3240                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3241      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3242 #endif
3243                 endif
3244               endif
3245               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3246             else
3247               ii=ii+1
3248               l_homo(k,ii)=.false.
3249             endif
3250             enddo
3251           enddo
3252         lim_odl=ii
3253         endif
3254 c
3255 c     Theta, dihedral and SC retraints
3256 c
3257         if (waga_angle.gt.0.0d0) then
3258 c         open (ientin,file=tpl_k_sigma_dih,status='old')
3259 c         do irec=1,maxres-3 ! loop for reading sigma_dih
3260 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3261 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3262 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3263 c    &                            sigma_dih(k,i+nnt-1)
3264 c         enddo
3265 c1402   continue
3266 c         close (ientin)
3267           do i = nnt+3,nct
3268             if (idomain(k,i).eq.0) then
3269                sigma_dih(k,i)=0.0
3270                cycle
3271             endif
3272             dih(k,i)=phiref(i) ! right?
3273 c           read (ientin,*) sigma_dih(k,i) ! original variant
3274 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
3275 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3276 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3277 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
3278 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
3279
3280             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3281      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3282 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3283 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
3284 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3285 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3286 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
3287 c   Instead of res sim other local measure of b/b str reliability possible
3288             if (sigma_dih(k,i).ne.0)
3289      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3290 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3291           enddo
3292           lim_dih=nct-nnt-2
3293         endif
3294
3295         if (waga_theta.gt.0.0d0) then
3296 c         open (ientin,file=tpl_k_sigma_theta,status='old')
3297 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3298 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3299 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3300 c    &                              sigma_theta(k,i+nnt-1)
3301 c         enddo
3302 c1403   continue
3303 c         close (ientin)
3304
3305           do i = nnt+2,nct ! right? without parallel.
3306 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
3307 c         do i=ithet_start,ithet_end ! with FG parallel.
3308              if (idomain(k,i).eq.0) then
3309               sigma_theta(k,i)=0.0
3310               cycle
3311              endif
3312              thetatpl(k,i)=thetaref(i)
3313 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3314 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
3315 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3316 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
3317 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
3318              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3319      &                        rescore(k,i-2))/3.0
3320 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3321              if (sigma_theta(k,i).ne.0)
3322      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3323
3324 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3325 c                             rescore(k,i-2) !  right expression ?
3326 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3327           enddo
3328         endif
3329
3330         if (waga_d.gt.0.0d0) then
3331 c       open (ientin,file=tpl_k_sigma_d,status='old')
3332 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3333 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3334 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3335 c    &                          sigma_d(k,i+nnt-1)
3336 c         enddo
3337 c1404   continue
3338
3339           do i = nnt,nct ! right? without parallel.
3340 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
3341 c         do i=loc_start,loc_end ! with FG parallel.
3342                if (itype(i).eq.10) cycle
3343                if (idomain(k,i).eq.0 ) then
3344                   sigma_d(k,i)=0.0
3345                   cycle
3346                endif
3347                xxtpl(k,i)=xxref(i)
3348                yytpl(k,i)=yyref(i)
3349                zztpl(k,i)=zzref(i)
3350 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3351 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3352 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3353 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
3354                sigma_d(k,i)=rescore3(k,i) !  right expression ?
3355                if (sigma_d(k,i).ne.0)
3356      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3357
3358 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
3359 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3360 c              read (ientin,*) sigma_d(k,i) ! 1st variant
3361           enddo
3362         endif
3363       enddo
3364 c
3365 c remove distance restraints not used in any model from the list
3366 c shift data in all arrays
3367 c
3368       if (waga_dist.ne.0.0d0) then
3369         ii=0
3370         liiflag=.true.
3371         do i=nnt,nct-2
3372          do j=i+2,nct
3373           ii=ii+1
3374           if (ii_in_use(ii).eq.0.and.liiflag) then
3375             liiflag=.false.
3376             iistart=ii
3377           endif
3378           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3379      &                   .not.liiflag.and.ii.eq.lim_odl) then
3380              if (ii.eq.lim_odl) then
3381               iishift=ii-iistart+1
3382              else
3383               iishift=ii-iistart
3384              endif
3385              liiflag=.true.
3386              do ki=iistart,lim_odl-iishift
3387               ires_homo(ki)=ires_homo(ki+iishift)
3388               jres_homo(ki)=jres_homo(ki+iishift)
3389               ii_in_use(ki)=ii_in_use(ki+iishift)
3390               do k=1,constr_homology
3391                odl(k,ki)=odl(k,ki+iishift)
3392                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3393                l_homo(k,ki)=l_homo(k,ki+iishift)
3394               enddo
3395              enddo
3396              ii=ii-iishift
3397              lim_odl=lim_odl-iishift
3398           endif
3399          enddo
3400         enddo
3401       endif
3402
3403       endif ! .not. klapaucjusz
3404
3405       if (constr_homology.gt.0) call homology_partition
3406       if (constr_homology.gt.0) call init_int_table
3407 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3408 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3409 c
3410 c Print restraints
3411 c
3412       if (.not.out_template_restr) return
3413 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3414       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3415        write (iout,*) "Distance restraints from templates"
3416        do ii=1,lim_odl
3417        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3418      &  ii,ires_homo(ii),jres_homo(ii),
3419      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3420      &  ki=1,constr_homology)
3421        enddo
3422        write (iout,*) "Dihedral angle restraints from templates"
3423        do i=nnt+3,nct
3424         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3425      &      (rad2deg*dih(ki,i),
3426      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3427        enddo
3428        write (iout,*) "Virtual-bond angle restraints from templates"
3429        do i=nnt+2,nct
3430         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3431      &      (rad2deg*thetatpl(ki,i),
3432      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3433        enddo
3434        write (iout,*) "SC restraints from templates"
3435        do i=nnt,nct
3436         write(iout,'(i5,100(4f8.2,4x))') i,
3437      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3438      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3439        enddo
3440       endif
3441 c -----------------------------------------------------------------
3442       return
3443       end
3444 c----------------------------------------------------------------------
3445 #ifdef WINIFL
3446       subroutine flush(iu)
3447       return
3448       end
3449 #endif
3450 #ifdef AIX
3451       subroutine flush(iu)
3452       call flush_(iu)
3453       return
3454       end
3455 #endif
3456 c------------------------------------------------------------------------------
3457       subroutine copy_to_tmp(source)
3458       include "DIMENSIONS"
3459       include "COMMON.IOUNITS"
3460       character*(*) source
3461       character* 256 tmpfile
3462       integer ilen
3463       external ilen
3464       logical ex
3465       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3466       inquire(file=tmpfile,exist=ex)
3467       if (ex) then
3468         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3469      &   " to temporary directory..."
3470         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3471         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3472       endif
3473       return
3474       end
3475 c------------------------------------------------------------------------------
3476       subroutine move_from_tmp(source)
3477       include "DIMENSIONS"
3478       include "COMMON.IOUNITS"
3479       character*(*) source
3480       integer ilen
3481       external ilen
3482       write (*,*) "Moving ",source(:ilen(source)),
3483      & " from temporary directory to working directory"
3484       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3485       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3486       return
3487       end
3488 c------------------------------------------------------------------------------
3489       subroutine random_init(seed)
3490 C
3491 C Initialize random number generator
3492 C
3493       implicit real*8 (a-h,o-z)
3494       include 'DIMENSIONS'
3495 #ifdef MPI
3496       include 'mpif.h'
3497       logical OKRandom, prng_restart
3498       real*8  r1
3499       integer iseed_array(4)
3500 #endif
3501       include 'COMMON.IOUNITS'
3502       include 'COMMON.TIME1'
3503       include 'COMMON.THREAD'
3504       include 'COMMON.SBRIDGE'
3505       include 'COMMON.CONTROL'
3506       include 'COMMON.MCM'
3507       include 'COMMON.MAP'
3508       include 'COMMON.HEADER'
3509       include 'COMMON.CSA'
3510       include 'COMMON.CHAIN'
3511       include 'COMMON.MUCA'
3512       include 'COMMON.MD'
3513       include 'COMMON.FFIELD'
3514       include 'COMMON.SETUP'
3515       iseed=-dint(dabs(seed))
3516       if (iseed.eq.0) then
3517         write (iout,'(/80(1h*)/20x,a/80(1h*))')
3518      &    'Random seed undefined. The program will stop.'
3519         write (*,'(/80(1h*)/20x,a/80(1h*))')
3520      &    'Random seed undefined. The program will stop.'
3521 #ifdef MPI
3522         call mpi_finalize(mpi_comm_world,ierr)
3523 #endif
3524         stop 'Bad random seed.'
3525       endif
3526 #ifdef MPI
3527       if (fg_rank.eq.0) then
3528       seed=seed*(me+1)+1
3529 #ifdef AMD64
3530       if(me.eq.king)
3531      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3532       OKRandom = prng_restart(me,iseed)
3533 #else
3534       do i=1,4
3535        tmp=65536.0d0**(4-i)
3536        iseed_array(i) = dint(seed/tmp)
3537        seed=seed-iseed_array(i)*tmp
3538       enddo
3539       if(me.eq.king)
3540      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3541      &                 (iseed_array(i),i=1,4)
3542       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3543      &                 (iseed_array(i),i=1,4)
3544       OKRandom = prng_restart(me,iseed_array)
3545 #endif
3546       if (OKRandom) then
3547 c        r1 = prng_next(me)
3548         r1=ran_number(0.0D0,1.0D0)
3549         if(me.eq.king)
3550      &   write (iout,*) 'ran_num',r1
3551         if (r1.lt.0.0d0) OKRandom=.false.
3552       endif
3553       if (.not.OKRandom) then
3554         write (iout,*) 'PRNG IS NOT WORKING!!!'
3555         print *,'PRNG IS NOT WORKING!!!'
3556         if (me.eq.0) then
3557          call flush(iout)
3558          call mpi_abort(mpi_comm_world,error_msg,ierr)
3559          stop
3560         else
3561          write (iout,*) 'too many processors for parallel prng'
3562          write (*,*) 'too many processors for parallel prng'
3563          call flush(iout)
3564          stop
3565         endif
3566       endif
3567       endif
3568 #else
3569       call vrndst(iseed)
3570       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3571 #endif
3572       return
3573       end
3574 c----------------------------------------------------------------------
3575       subroutine read_klapaucjusz
3576
3577       include 'DIMENSIONS'
3578 #ifdef MPI
3579       include 'mpif.h'
3580 #endif
3581       include 'COMMON.SETUP'
3582       include 'COMMON.CONTROL'
3583       include 'COMMON.CHAIN'
3584       include 'COMMON.IOUNITS'
3585       include 'COMMON.MD'
3586       include 'COMMON.GEO'
3587       include 'COMMON.INTERACT'
3588       include 'COMMON.NAMES'
3589       character*256 fragfile
3590       integer ninclust(maxclust),inclust(max_template,maxclust),
3591      &  nresclust(maxclust),iresclust(maxres,maxclust)
3592
3593       character*2 kic2
3594       character*24 model_ki_dist, model_ki_angle
3595       character*500 controlcard
3596       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
3597       logical lprn /.true./
3598       integer ilen
3599       external ilen
3600       logical liiflag
3601 c
3602 c
3603       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3604       double precision, dimension (max_template,maxres) :: rescore
3605       double precision, dimension (max_template,maxres) :: rescore2
3606       character*24 pdbfile,tpl_k_rescore
3607
3608 c
3609 c For new homol impl
3610 c
3611       include 'COMMON.VAR'
3612 c
3613       call getenv("FRAGFILE",fragfile)
3614       open(ientin,file=fragfile,status="old",err=10)
3615       read(ientin,*) constr_homology,nclust
3616       l_homo = .false.
3617       sigma_theta=0.0
3618       sigma_d=0.0
3619       sigma_dih=0.0
3620 c Read pdb files
3621       do k=1,constr_homology
3622         read(ientin,'(a)') pdbfile
3623         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3624      &  pdbfile(:ilen(pdbfile))
3625         open(ipdbin,file=pdbfile,status='old',err=33)
3626         goto 34
3627   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3628      &  pdbfile(:ilen(pdbfile))
3629         stop
3630   34    continue
3631         unres_pdb=.false.
3632         call readpdb_template(k)
3633         do i=1,nres
3634           rescore(k,i)=0.2d0
3635           rescore2(k,i)=1.0d0
3636         enddo
3637       enddo
3638 c Read clusters
3639       do i=1,nclust
3640         read(ientin,*) ninclust(i),nresclust(i)
3641         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3642         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3643       enddo
3644 c
3645 c Loop over clusters
3646 c
3647       do l=1,nclust
3648         do ll = 1,ninclust(l)
3649
3650         k = inclust(ll,l)
3651         do i=1,nres
3652           idomain(k,i)=0
3653         enddo
3654         do i=1,nresclust(l)
3655           if (nnt.gt.1)  then
3656             idomain(k,iresclust(i,l)+1) = 1
3657           else
3658             idomain(k,iresclust(i,l)) = 1
3659           endif
3660         enddo
3661 c
3662 c     Distance restraints
3663 c
3664 c          ... --> odl(k,ii)
3665 C Copy the coordinates from reference coordinates (?)
3666         do i=1,2*nres
3667           do j=1,3
3668             c(j,i)=chomo(j,i,k)
3669 c           write (iout,*) "c(",j,i,") =",c(j,i)
3670           enddo
3671         enddo
3672         call int_from_cart(.true.,.false.)
3673         call sc_loc_geom(.true.)
3674         do i=1,nres
3675           thetaref(i)=theta(i)
3676           phiref(i)=phi(i)
3677         enddo
3678         if (waga_dist.ne.0.0d0) then
3679           ii=0
3680           do i = nnt,nct-2
3681             do j=i+2,nct
3682
3683               x12=c(1,i)-c(1,j)
3684               y12=c(2,i)-c(2,j)
3685               z12=c(3,i)-c(3,j)
3686               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3687 c              write (iout,*) k,i,j,distal,dist2_cut
3688
3689             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3690      &            .and. distal.le.dist2_cut ) then
3691
3692               ii=ii+1
3693               ii_in_use(ii)=1
3694               l_homo(k,ii)=.true.
3695
3696 c             write (iout,*) "k",k
3697 c             write (iout,*) "i",i," j",j," constr_homology",
3698 c    &                       constr_homology
3699               ires_homo(ii)=i
3700               jres_homo(ii)=j
3701               odl(k,ii)=distal
3702               if (read2sigma) then
3703                 sigma_odl(k,ii)=0
3704                 do ik=i,j
3705                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3706                 enddo
3707                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3708                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3709      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3710               else
3711                 if (odl(k,ii).le.dist_cut) then
3712                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3713                 else
3714 #ifdef OLDSIGMA
3715                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3716      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3717 #else
3718                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3719      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3720 #endif
3721                 endif
3722               endif
3723               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3724             else
3725               ii=ii+1
3726 c              l_homo(k,ii)=.false.
3727             endif
3728             enddo
3729           enddo
3730         lim_odl=ii
3731         endif
3732 c
3733 c     Theta, dihedral and SC retraints
3734 c
3735         if (waga_angle.gt.0.0d0) then
3736           do i = nnt+3,nct
3737             if (idomain(k,i).eq.0) then
3738 c               sigma_dih(k,i)=0.0
3739                cycle
3740             endif
3741             dih(k,i)=phiref(i)
3742             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3743      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3744 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3745 c     &       " sigma_dihed",sigma_dih(k,i)
3746             if (sigma_dih(k,i).ne.0)
3747      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3748           enddo
3749           lim_dih=nct-nnt-2
3750         endif
3751
3752         if (waga_theta.gt.0.0d0) then
3753           do i = nnt+2,nct
3754              if (idomain(k,i).eq.0) then
3755 c              sigma_theta(k,i)=0.0
3756               cycle
3757              endif
3758              thetatpl(k,i)=thetaref(i)
3759              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3760      &                        rescore(k,i-2))/3.0
3761              if (sigma_theta(k,i).ne.0)
3762      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3763           enddo
3764         endif
3765
3766         if (waga_d.gt.0.0d0) then
3767           do i = nnt,nct
3768                if (itype(i).eq.10) cycle
3769                if (idomain(k,i).eq.0 ) then
3770 c                  sigma_d(k,i)=0.0
3771                   cycle
3772                endif
3773                xxtpl(k,i)=xxref(i)
3774                yytpl(k,i)=yyref(i)
3775                zztpl(k,i)=zzref(i)
3776                sigma_d(k,i)=rescore(k,i)
3777                if (sigma_d(k,i).ne.0)
3778      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3779                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3780           enddo
3781         endif
3782       enddo ! l
3783       enddo ! ll
3784 c
3785 c remove distance restraints not used in any model from the list
3786 c shift data in all arrays
3787 c
3788       if (waga_dist.ne.0.0d0) then
3789         ii=0
3790         liiflag=.true.
3791         do i=nnt,nct-2
3792          do j=i+2,nct
3793           ii=ii+1
3794           if (ii_in_use(ii).eq.0.and.liiflag) then
3795             liiflag=.false.
3796             iistart=ii
3797           endif
3798           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3799      &                   .not.liiflag.and.ii.eq.lim_odl) then
3800              if (ii.eq.lim_odl) then
3801               iishift=ii-iistart+1
3802              else
3803               iishift=ii-iistart
3804              endif
3805              liiflag=.true.
3806              do ki=iistart,lim_odl-iishift
3807               ires_homo(ki)=ires_homo(ki+iishift)
3808               jres_homo(ki)=jres_homo(ki+iishift)
3809               ii_in_use(ki)=ii_in_use(ki+iishift)
3810               do k=1,constr_homology
3811                odl(k,ki)=odl(k,ki+iishift)
3812                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3813                l_homo(k,ki)=l_homo(k,ki+iishift)
3814               enddo
3815              enddo
3816              ii=ii-iishift
3817              lim_odl=lim_odl-iishift
3818           endif
3819          enddo
3820         enddo
3821       endif
3822
3823       return
3824    10 stop "Error in fragment file"
3825       end