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