source/unres/src-HCD-5D/readrtns_CSA.F

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