8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
14 C Read job setup parameters
16 C Read force-field parameters except weights
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
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
37 cfmc if (modecalc.eq.10) call mcmfread
38 c Read energy-term weights and disulfide parameters
40 C Read molecule information, molecule geometry, energy-term weights, and
41 C restraints if requested
43 C Print restraint information
45 if (.not. out1file .or. me.eq.king) 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",
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),
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"
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)
72 write (iout,*) "The ipeak array"
74 write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
80 c print *,"Processor",myrank," leaves READRTNS"
83 C-------------------------------------------------------------------------------
84 subroutine read_control
92 logical OKRandom, prng_restart
95 include 'COMMON.IOUNITS'
96 include 'COMMON.TIME1'
97 include 'COMMON.THREAD'
98 include 'COMMON.SBRIDGE'
99 include 'COMMON.CONTROL'
100 include 'COMMON.SAXS'
103 include 'COMMON.HEADER'
105 include 'COMMON.CHAIN'
106 include 'COMMON.MUCA'
108 include 'COMMON.FFIELD'
109 include 'COMMON.INTERACT'
110 include 'COMMON.SETUP'
111 include 'COMMON.SPLITELE'
112 include 'COMMON.SHIELD'
115 integer KDIAG,ICORFL,IXDR
116 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
117 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
119 character*320 controlcard
120 double precision seed
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
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)
172 safety = 60.0d0*safety
174 call reada(controlcard,"T_BATH",t_bath,300.0d0)
175 minim=(index(controlcard,'MINIMIZE').gt.0)
176 dccart=(index(controlcard,'CART').gt.0)
177 overlapsc=(index(controlcard,'OVERLAP').gt.0)
178 overlapsc=.not.overlapsc
179 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
180 searchsc=.not.searchsc
181 sideadd=(index(controlcard,'SIDEADD').gt.0)
182 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
183 mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
184 outpdb=(index(controlcard,'PDBOUT').gt.0)
185 outmol2=(index(controlcard,'MOL2OUT').gt.0)
186 pdbref=(index(controlcard,'PDBREF').gt.0)
187 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
188 indpdb=index(controlcard,'PDBSTART')
189 extconf=(index(controlcard,'EXTCONF').gt.0)
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
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)
221 if (index(controlcard,"CASC").gt.0) then
223 c else if (index(controlcard,"CAONLY").gt.0) then
225 else if (index(controlcard,"SCONLY").gt.0) then
229 c write (iout,*) "ERROR! Must specify CASC CAONLY or SCONLY when",
230 c & " specifying REFSTR, PDBREF or REGULAR."
234 if (index(controlcard,'CHECKGRAD').gt.0) then
236 if (index(controlcard,' CART').gt.0) then
238 elseif (index(controlcard,'CARINT').gt.0) then
243 call reada(controlcard,'DELTA',aincr,1.0d-7)
244 c write (iout,*) "icheckgrad",icheckgrad
245 elseif (index(controlcard,'THREAD').gt.0) then
247 call readi(controlcard,'THREAD',nthread,0)
248 if (nthread.gt.0) then
249 call reada(controlcard,'WEIDIS',weidis,0.1D0)
252 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
253 stop 'Error termination in Read_Control.'
255 else if (index(controlcard,'MCMA').gt.0) then
257 else if (index(controlcard,'MCEE').gt.0) then
259 else if (index(controlcard,'MULTCONF').gt.0) then
261 else if (index(controlcard,'MAP').gt.0) then
263 call readi(controlcard,'MAP',nmap,0)
264 else if (index(controlcard,'CSA').gt.0) then
266 crc else if (index(controlcard,'ZSCORE').gt.0) then
268 crc ZSCORE is rm from UNRES, modecalc=9 is available
271 cfcm else if (index(controlcard,'MCMF').gt.0) then
273 else if (index(controlcard,'SOFTREG').gt.0) then
275 else if (index(controlcard,'CHECK_BOND').gt.0) then
277 else if (index(controlcard,'TEST').gt.0) then
279 else if (index(controlcard,'MD').gt.0) then
281 else if (index(controlcard,'RE ').gt.0) then
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 ))
290 write (iout,*) 'MUCADYN=',mucadyn
291 write (iout,*) 'MUCASMOOTH=',muca_smooth
294 iscode=index(controlcard,'ONE_LETTER')
295 indphi=index(controlcard,'PHI')
296 indback=index(controlcard,'BACK')
297 iranconf=index(controlcard,'RAND_CONF')
298 i2ndstr=index(controlcard,'USE_SEC_PRED')
299 gradout=index(controlcard,'GRADOUT').gt.0
300 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
301 C DISTCHAINMAX become obsolete for periodic boundry condition
302 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
303 C Reading the dimensions of box in x,y,z coordinates
304 call reada(controlcard,'BOXX',boxxsize,100.0d0)
305 call reada(controlcard,'BOXY',boxysize,100.0d0)
306 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
307 c Cutoff range for interactions
308 call reada(controlcard,"R_CUT",r_cut,15.0d0)
309 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
310 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
311 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
312 if (lipthick.gt.0.0d0) then
313 bordliptop=(boxzsize+lipthick)/2.0
314 bordlipbot=bordliptop-lipthick
316 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
317 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
318 buflipbot=bordlipbot+lipbufthick
319 bufliptop=bordliptop-lipbufthick
320 if ((lipbufthick*2.0d0).gt.lipthick)
321 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
323 write(iout,*) "bordliptop=",bordliptop
324 write(iout,*) "bordlipbot=",bordlipbot
325 write(iout,*) "bufliptop=",bufliptop
326 write(iout,*) "buflipbot=",buflipbot
327 write (iout,*) "SHIELD MODE",shield_mode
328 if (TUBElog.gt.0) then
329 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
330 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
331 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
332 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
333 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
334 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
335 buftubebot=bordtubebot+tubebufthick
336 buftubetop=bordtubetop-tubebufthick
338 if (me.eq.king .or. .not.out1file )
339 & write (iout,*) "DISTCHAINMAX",distchainmax
341 if(me.eq.king.or..not.out1file)
342 & write (iout,'(2a)') diagmeth(kdiag),
343 & ' routine used to diagonalize matrices.'
346 c--------------------------------------------------------------------------
347 subroutine read_REMDpar
353 include 'COMMON.IOUNITS'
354 include 'COMMON.TIME1'
357 include 'COMMON.LANGEVIN'
360 include 'COMMON.LANGEVIN.lang0.5diag'
362 include 'COMMON.LANGEVIN.lang0'
365 include 'COMMON.INTERACT'
366 include 'COMMON.NAMES'
368 include 'COMMON.REMD'
369 include 'COMMON.CONTROL'
370 include 'COMMON.SETUP'
372 character*320 controlcard
373 character*3200 controlcard1
374 integer iremd_m_total,i
376 if(me.eq.king.or..not.out1file)
377 & write (iout,*) "REMD setup"
379 call card_concat(controlcard)
380 call readi(controlcard,"NREP",nrep,3)
381 call readi(controlcard,"NSTEX",nstex,1000)
382 call reada(controlcard,"RETMIN",retmin,10.0d0)
383 call reada(controlcard,"RETMAX",retmax,1000.0d0)
384 mremdsync=(index(controlcard,'SYNC').gt.0)
385 call readi(controlcard,"NSYN",i_sync_step,100)
386 restart1file=(index(controlcard,'REST1FILE').gt.0)
387 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
388 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
389 if(max_cache_traj_use.gt.max_cache_traj)
390 & max_cache_traj_use=max_cache_traj
391 if(me.eq.king.or..not.out1file) then
392 cd if (traj1file) then
393 crc caching is in testing - NTWX is not ignored
394 cd write (iout,*) "NTWX value is ignored"
395 cd write (iout,*) " trajectory is stored to one file by master"
396 cd write (iout,*) " before exchange at NSTEX intervals"
398 write (iout,*) "NREP= ",nrep
399 write (iout,*) "NSTEX= ",nstex
400 write (iout,*) "SYNC= ",mremdsync
401 write (iout,*) "NSYN= ",i_sync_step
402 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
405 if (index(controlcard,'TLIST').gt.0) then
407 call card_concat(controlcard1)
408 read(controlcard1,*) (remd_t(i),i=1,nrep)
409 if(me.eq.king.or..not.out1file)
410 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
413 if (index(controlcard,'MLIST').gt.0) then
415 call card_concat(controlcard1)
416 read(controlcard1,*) (remd_m(i),i=1,nrep)
417 if(me.eq.king.or..not.out1file) then
418 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
421 iremd_m_total=iremd_m_total+remd_m(i)
423 write (iout,*) 'Total number of replicas ',iremd_m_total
428 & "Adaptive (PMF-biased) umbrella sampling will be run"
431 if(me.eq.king.or..not.out1file)
432 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
435 c--------------------------------------------------------------------------
436 subroutine read_MDpar
442 include 'COMMON.IOUNITS'
443 include 'COMMON.TIME1'
445 include 'COMMON.QRESTR'
447 include 'COMMON.LANGEVIN'
450 include 'COMMON.LANGEVIN.lang0.5diag'
452 include 'COMMON.LANGEVIN.lang0'
455 include 'COMMON.INTERACT'
456 include 'COMMON.NAMES'
458 include 'COMMON.SETUP'
459 include 'COMMON.CONTROL'
460 include 'COMMON.SPLITELE'
461 include 'COMMON.FFIELD'
463 character*320 controlcard
467 call card_concat(controlcard)
468 call readi(controlcard,"NSTEP",n_timestep,1000000)
469 call readi(controlcard,"NTWE",ntwe,100)
470 call readi(controlcard,"NTWX",ntwx,1000)
471 call reada(controlcard,"DT",d_time,1.0d-1)
472 call reada(controlcard,"DVMAX",dvmax,2.0d1)
473 call reada(controlcard,"DAMAX",damax,1.0d1)
474 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
475 call readi(controlcard,"LANG",lang,0)
476 RESPA = index(controlcard,"RESPA") .gt. 0
477 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
478 ntime_split0=ntime_split
479 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
480 ntime_split0=ntime_split
481 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
482 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
483 rest = index(controlcard,"REST").gt.0
484 tbf = index(controlcard,"TBF").gt.0
485 usampl = index(controlcard,"USAMPL").gt.0
486 scale_umb = index(controlcard,"SCALE_UMB").gt.0
487 adaptive = index(controlcard,"ADAPTIVE").gt.0
488 mdpdb = index(controlcard,"MDPDB").gt.0
489 call reada(controlcard,"T_BATH",t_bath,300.0d0)
490 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
491 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
492 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
493 if (count_reset_moment.eq.0) count_reset_moment=1000000000
494 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
495 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
496 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
497 if (count_reset_vel.eq.0) count_reset_vel=1000000000
498 large = index(controlcard,"LARGE").gt.0
499 print_compon = index(controlcard,"PRINT_COMPON").gt.0
500 rattle = index(controlcard,"RATTLE").gt.0
501 preminim = index(controlcard,"PREMINIM").gt.0
502 if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then
506 write (iout,*) "PREMINIM ",preminim
508 if (index(controlcard,'CART').gt.0) dccart=.true.
509 write (iout,*) "dccart ",dccart
510 write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR')
511 if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim
513 c if performing umbrella sampling, fragments constrained are read from the fragment file
516 write (iout,*) "Umbrella sampling will be run"
517 if (scale_umb.and.adaptive) then
518 write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
519 write (iout,*) "Select one of those and re-run the job."
522 if (scale_umb) write (iout,*)
523 &"Umbrella-restraint force constants will be scaled by temperature"
527 if(me.eq.king.or..not.out1file) then
529 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
531 write (iout,'(a)') "The units are:"
532 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
533 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
534 & " acceleration: angstrom/(48.9 fs)**2"
535 write (iout,'(a)') "energy: kcal/mol, temperature: K"
537 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
538 write (iout,'(a60,f10.5,a)')
539 & "Initial time step of numerical integration:",d_time,
541 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
543 write (iout,'(2a,i4,a)')
544 & "A-MTS algorithm used; initial time step for fast-varying",
545 & " short-range forces split into",ntime_split," steps."
546 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
547 & r_cut," lambda",rlamb
549 write (iout,'(2a,f10.5)')
550 & "Maximum acceleration threshold to reduce the time step",
551 & "/increase split number:",damax
552 write (iout,'(2a,f10.5)')
553 & "Maximum predicted energy drift to reduce the timestep",
554 & "/increase split number:",edriftmax
555 write (iout,'(a60,f10.5)')
556 & "Maximum velocity threshold to reduce velocities:",dvmax
557 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
558 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
559 if (rattle) write (iout,'(a60)')
560 & "Rattle algorithm used to constrain the virtual bonds"
564 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
565 call reada(controlcard,"RWAT",rwat,1.4d0)
566 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
567 surfarea=index(controlcard,"SURFAREA").gt.0
568 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
569 if(me.eq.king.or..not.out1file)then
570 write (iout,'(/a,$)') "Langevin dynamics calculation"
573 & " with direct integration of Langevin equations"
574 else if (lang.eq.2) then
575 write (iout,'(a/)') " with TINKER stochasic MD integrator"
576 else if (lang.eq.3) then
577 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
578 else if (lang.eq.4) then
579 write (iout,'(a/)') " in overdamped mode"
581 write (iout,'(//a,i5)')
582 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
585 write (iout,'(a60,f10.5)') "Temperature:",t_bath
586 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
587 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
588 write (iout,'(a60,f10.5)')
589 & "Scaling factor of the friction forces:",scal_fric
590 if (surfarea) write (iout,'(2a,i10,a)')
591 & "Friction coefficients will be scaled by solvent-accessible",
592 & " surface area every",reset_fricmat," steps."
594 c Calculate friction coefficients and bounds of stochastic forces
595 eta=6*pi*cPoise*etawat
596 if(me.eq.king.or..not.out1file)
597 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
599 gamp=scal_fric*(pstok+rwat)*eta
600 stdfp=dsqrt(2*Rb*t_bath/d_time)
602 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
603 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
605 if(me.eq.king.or..not.out1file)then
606 write (iout,'(/2a/)')
607 & "Radii of site types and friction coefficients and std's of",
608 & " stochastic forces of fully exposed sites"
609 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
611 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
612 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
616 if(me.eq.king.or..not.out1file)then
617 write (iout,'(a)') "Berendsen bath calculation"
618 write (iout,'(a60,f10.5)') "Temperature:",t_bath
619 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
621 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
622 & count_reset_moment," steps"
624 & write (iout,'(a,i10,a)')
625 & "Velocities will be reset at random every",count_reset_vel,
629 if(me.eq.king.or..not.out1file)
630 & write (iout,'(a31)') "Microcanonical mode calculation"
632 if(me.eq.king.or..not.out1file)then
633 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
635 write(iout,*) "MD running with constraints."
636 write(iout,*) "Equilibration time ", eq_time, " mtus."
637 write(iout,*) "Constraining ", nfrag," fragments."
638 write(iout,*) "Length of each fragment, weight and q0:"
640 write (iout,*) "Set of restraints #",iset
642 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
643 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
645 write(iout,*) "constraints between ", npair, "fragments."
646 write(iout,*) "constraint pairs, weights and q0:"
648 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
649 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
651 write(iout,*) "angle constraints within ", nfrag_back,
652 & "backbone fragments."
654 write(iout,*) "fragment, weights, q0:"
656 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
657 & ifrag_back(2,i,iset),
658 & wfrag_back(1,i,iset),qin_back(1,i,iset),
659 & wfrag_back(2,i,iset),qin_back(2,i,iset),
660 & wfrag_back(3,i,iset),qin_back(3,i,iset)
663 write(iout,*) "fragment, weights:"
665 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
666 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
667 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
671 iset=mod(kolor,nset)+1
674 if(me.eq.king.or..not.out1file)
675 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
678 c------------------------------------------------------------------------------
681 C Read molecular data.
687 integer error_msg,ierror,ierr,ierrcode
689 include 'COMMON.IOUNITS'
692 include 'COMMON.INTERACT'
693 include 'COMMON.LOCAL'
694 include 'COMMON.NAMES'
695 include 'COMMON.CHAIN'
696 include 'COMMON.FFIELD'
697 include 'COMMON.SBRIDGE'
698 include 'COMMON.HEADER'
699 include 'COMMON.CONTROL'
700 include 'COMMON.SAXS'
701 include 'COMMON.DBASE'
702 include 'COMMON.THREAD'
703 include 'COMMON.CONTACTS'
704 include 'COMMON.TORCNSTR'
705 include 'COMMON.TIME1'
706 include 'COMMON.BOUNDS'
708 include 'COMMON.SETUP'
709 include 'COMMON.SHIELD'
710 character*4 sequence(maxres)
712 double precision x(maxvar)
713 character*256 pdbfile
714 character*400 weightcard
715 character*80 weightcard_t,ucase
716 integer itype_pdb(maxres)
717 common /pizda/ itype_pdb
718 logical seq_comp,fail
719 double precision energia(0:n_ene)
720 double precision secprob(3,maxdih_constr)
722 double precision phihel,phibet,sigmahel,sigmabet
723 integer iti,nsi,maxsi
727 integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2
728 double precision sumv
730 C Read PDB structure if applicable
732 if (indpdb.gt.0 .or. pdbref) then
733 read(inp,'(a)') pdbfile
734 if(me.eq.king.or..not.out1file)
735 & write (iout,'(2a)') 'PDB data will be read from file ',
736 & pdbfile(:ilen(pdbfile))
737 open(ipdbin,file=pdbfile,status='old',err=33)
739 33 write (iout,'(a)') 'Error opening PDB file.'
750 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
751 & (crefjlee(j,i+nres),j=1,3)
754 if(me.eq.king.or..not.out1file)
755 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
756 & ' nstart_sup=',nstart_sup
758 itype_pdb(i)=itype(i)
762 nct=nstart_sup+nsup-1
763 call contact(.false.,ncont_ref,icont_ref,co)
766 if(me.eq.king.or..not.out1file)
767 & write(iout,*)'Adding sidechains'
771 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
774 do while (fail.and.nsi.le.maxsi)
775 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
778 if(fail) write(iout,*)'Adding sidechain failed for res ',
779 & i,' after ',nsi,' trials'
784 if (indpdb.eq.0) then
785 C Read sequence if not taken from the pdb file.
787 c print *,'nres=',nres
788 if (iscode.gt.0) then
789 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
791 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
793 C Convert sequence to numeric code
795 itype(i)=rescode(i,sequence(i),iscode)
799 c print '(20i4)',(itype(i),i=1,nres)
802 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
804 if (itype(i).eq.ntyp1) then
808 else if (iabs(itype(i+1)).ne.20) then
810 else if (iabs(itype(i)).ne.20) then
817 if(me.eq.king.or..not.out1file)then
818 write (iout,*) "ITEL"
820 write (iout,*) i,itype(i),itel(i)
822 print *,'Call Read_Bridge.'
826 cd print *,'NNT=',NNT,' NCT=',NCT
827 call seq2chains(nres,itype,nchain,chain_length,chain_border,
830 chain_border1(2,1)=chain_border(2,1)+1
832 chain_border1(1,i)=chain_border(1,i)-1
833 chain_border1(2,i)=chain_border(2,i)+1
835 chain_border1(1,nchain)=chain_border(1,nchain)-1
836 chain_border1(2,nchain)=nres
837 write(iout,*) "nres",nres," nchain",nchain
839 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
840 & chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
842 call chain_symmetry(nchain,nres,itype,chain_border,
843 & chain_length,npermchain,tabpermchain)
845 c write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
848 write(iout,*) "residue permutations"
850 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
853 if (itype(1).eq.ntyp1) nnt=2
854 if (itype(nres).eq.ntyp1) nct=nct-1
855 write (iout,*) "nnt",nnt," nct",nct
858 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
859 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
861 print*, 'init_dfa_vars finished!'
863 print*, 'read_dfa_info finished!'
867 if(me.eq.king.or..not.out1file)
868 & write (iout,'(a,i3)') 'nsup=',nsup
870 if (nsup.le.(nct-nnt+1)) then
871 do i=0,nct-nnt+1-nsup
872 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
878 & 'Error - sequences to be superposed do not match.'
881 do i=0,nsup-(nct-nnt+1)
882 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
884 nstart_sup=nstart_sup+i
890 & 'Error - sequences to be superposed do not match.'
893 if (nsup.eq.0) nsup=nct-nnt
894 if (nstart_sup.eq.0) nstart_sup=nnt
895 if (nstart_seq.eq.0) nstart_seq=nnt
896 if(me.eq.king.or..not.out1file)
897 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
898 & ' nstart_seq=',nstart_seq
901 C 8/13/98 Set limits to generating the dihedral angles
906 read (inp,*) ndih_constr
907 write (iout,*) "ndih_constr",ndih_constr
908 if (ndih_constr.gt.0) then
911 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
914 if(me.eq.king.or..not.out1file)then
917 & 'There are',ndih_constr,' restraints on gamma angles.'
919 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
926 phi0(i)=deg2rad*phi0(i)
927 drange(i)=deg2rad*drange(i)
929 C if(me.eq.king.or..not.out1file)
930 C & write (iout,*) 'FTORS',ftors
933 phibound(1,ii) = phi0(i)-drange(i)
934 phibound(2,ii) = phi0(i)+drange(i)
937 if (me.eq.king .or. .not.out1file) then
939 write (iout,'(a)') 'Boundaries in gamma angle sampling:'
941 write (iout,'(a3,i5,2f10.1)')
942 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
947 else if (ndih_constr.lt.0) then
949 call card_concat(weightcard)
950 call reada(weightcard,"PHIHEL",phihel,50.0D0)
951 call reada(weightcard,"PHIBET",phibet,180.0D0)
952 call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
953 call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
954 call reada(weightcard,"WDIHC",wdihc,0.591D0)
955 write (iout,*) "Weight of dihedral angle restraints",wdihc
956 read(inp,'(9x,3f7.3)')
957 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
958 write (iout,*) "The secprob array"
960 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
964 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
965 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
966 ndih_constr=ndih_constr+1
967 idih_constr(ndih_constr)=i
970 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
971 sumv=sumv+vpsipred(j,ndih_constr)
974 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
976 phibound(1,ndih_constr)=phihel*deg2rad
977 phibound(2,ndih_constr)=phibet*deg2rad
978 sdihed(1,ndih_constr)=sigmahel*deg2rad
979 sdihed(2,ndih_constr)=sigmabet*deg2rad
983 if(me.eq.king.or..not.out1file)then
986 & 'There are',ndih_constr,
987 & ' bimodal restraints on gamma angles.'
989 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
990 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
991 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
992 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
993 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
994 & (vpsipred(j,i),j=1,3)
1001 C first setting the theta boundaries to 0 to pi
1002 C this mean that there is no energy penalty for any angle occuring this can be applied
1003 C for generate random conformation but is not implemented in this way
1006 C thetabound(2,i)=pi
1008 C begin reading theta constrains this is quartic constrains allowing to
1009 C have smooth second derivative
1010 if (with_theta_constr) then
1011 C with_theta_constr is keyword allowing for occurance of theta constrains
1012 read (inp,*) ntheta_constr
1013 C ntheta_constr is the number of theta constrains
1014 if (ntheta_constr.gt.0) then
1015 C read (inp,*) ftors
1016 read (inp,*) (itheta_constr(i),theta_constr0(i),
1017 & theta_drange(i),for_thet_constr(i),
1018 & i=1,ntheta_constr)
1019 C the above code reads from 1 to ntheta_constr
1020 C itheta_constr(i) residue i for which is theta_constr
1021 C theta_constr0 the global minimum value
1022 C theta_drange is range for which there is no energy penalty
1023 C for_thet_constr is the force constant for quartic energy penalty
1025 if(me.eq.king.or..not.out1file)then
1027 & 'There are',ntheta_constr,' constraints on phi angles.'
1028 do i=1,ntheta_constr
1029 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1031 & for_thet_constr(i)
1034 do i=1,ntheta_constr
1035 theta_constr0(i)=deg2rad*theta_constr0(i)
1036 theta_drange(i)=deg2rad*theta_drange(i)
1038 C if(me.eq.king.or..not.out1file)
1039 C & write (iout,*) 'FTORS',ftors
1040 C do i=1,ntheta_constr
1041 C ii = itheta_constr(i)
1042 C thetabound(1,ii) = phi0(i)-drange(i)
1043 C thetabound(2,ii) = phi0(i)+drange(i)
1045 endif ! ntheta_constr.gt.0
1046 endif! with_theta_constr
1048 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1049 C write (iout,*) "with_dihed_constr ",with_dihed_constr
1050 c--- Zscore rms -------
1051 if (nz_start.eq.0) nz_start=nnt
1052 if (nz_end.eq.0 .and. nsup.gt.0) then
1054 else if (nz_end.eq.0) then
1057 if(me.eq.king.or..not.out1file)then
1058 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1059 write (iout,*) 'IZ_SC=',iz_sc
1061 c----------------------
1064 if (.not.pdbref) then
1065 call read_angles(inp,*38)
1068 38 write (iout,'(a)') 'Error reading reference structure.'
1070 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1071 stop 'Error reading reference structure'
1073 39 call chainbuild_extconf
1075 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1084 call contact(.true.,ncont_ref,icont_ref,co)
1088 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1090 if (constr_dist.gt.0) call read_dist_constr
1091 write (iout,*) "After read_dist_constr nhpb",nhpb
1092 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1094 call NMRpeak_partition
1095 if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1097 if(me.eq.king.or..not.out1file)
1098 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1101 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1103 if(me.eq.king.or..not.out1file)
1104 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1105 & icont_ref(1,i),' ',
1106 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1109 write (iout,*) "calling read_saxs_consrtr",nsaxs
1110 if (nsaxs.gt.0) call read_saxs_constr
1112 if (constr_homology.gt.0) then
1113 call read_constr_homology
1114 if (indpdb.gt.0 .or. pdbref) then
1117 c(j,i)=crefjlee(j,i)
1118 cref(j,i)=crefjlee(j,i)
1123 write (iout,*) "sc_loc_geom: Array C"
1125 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1126 & (c(j,i+nres),j=1,3)
1128 write (iout,*) "Array Cref"
1130 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1131 & (cref(j,i+nres),j=1,3)
1134 call int_from_cart1(.false.)
1135 call sc_loc_geom(.false.)
1137 thetaref(i)=theta(i)
1142 dc(j,i)=c(j,i+1)-c(j,i)
1143 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1148 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1149 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1158 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1159 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1160 & modecalc.ne.10) then
1161 C If input structure hasn't been supplied from the PDB file read or generate
1163 if (iranconf.eq.0 .and. .not. extconf) then
1164 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1165 & write (iout,'(a)') 'Initial geometry will be read in.'
1167 read(inp,'(8f10.5)',end=36,err=36)
1168 & ((c(l,k),l=1,3),k=1,nres),
1169 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1170 write (iout,*) "Exit READ_CART"
1171 c write (iout,'(8f10.5)')
1172 c & ((c(l,k),l=1,3),k=1,nres),
1173 c & ((c(l,k+nres),l=1,3),k=nnt,nct)
1177 dc(j,i)=c(j,i+1)-c(j,i)
1178 c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1182 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1184 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1185 c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1190 c dc_norm(j,i+nres)=0.0d0
1194 call int_from_cart1(.true.)
1195 write (iout,*) "Finish INT_TO_CART"
1196 c write (iout,*) "DC"
1198 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1199 c & (dc(j,i+nres),j=1,3)
1205 call read_angles(inp,*36)
1207 call chainbuild_extconf
1210 36 write (iout,'(a)') 'Error reading angle file.'
1212 call mpi_finalize( MPI_COMM_WORLD,IERR )
1214 stop 'Error reading angle file.'
1216 else if (extconf) then
1217 if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1218 & write (iout,'(a)') 'Extended chain initial geometry.'
1220 theta(i)=90d0*deg2rad
1223 phi(i)=180d0*deg2rad
1226 alph(i)=110d0*deg2rad
1229 omeg(i)=-120d0*deg2rad
1230 if (itype(i).le.0) omeg(i)=-omeg(i)
1233 call chainbuild_extconf
1235 if(me.eq.king.or..not.out1file)
1236 & write (iout,'(a)') 'Random-generated initial geometry.'
1239 if (me.eq.king .or. fg_rank.eq.0 .and. (
1240 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1244 call gen_rand_conf(itmp,*30)
1246 30 write (iout,*) 'Failed to generate random conformation',
1247 & ', itrial=',itrial
1248 write (*,*) 'Processor:',me,
1249 & ' Failed to generate random conformation',
1259 write (iout,'(a,i3,a)') 'Processor:',me,
1260 & ' error in generating random conformation.'
1261 write (*,'(a,i3,a)') 'Processor:',me,
1262 & ' error in generating random conformation.'
1265 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1270 & ' error in generating random conformation.'
1275 elseif (modecalc.eq.4) then
1276 read (inp,'(a)') intinname
1277 open (intin,file=intinname,status='old',err=333)
1278 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1279 & write (iout,'(a)') 'intinname',intinname
1280 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1282 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1284 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1286 stop 'Error opening angle file.'
1290 C Generate distance constraints, if the PDB structure is to be regularized.
1291 if (nthread.gt.0) then
1292 call read_threadbase
1295 if (me.eq.king .or. .not. out1file)
1297 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1298 write (iout,'(/a,i3,a)')
1299 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1300 write (iout,'(20i4)') (iss(i),i=1,ns)
1302 write(iout,*)"Running with dynamic disulfide-bond formation"
1304 write (iout,'(/a/)') 'Pre-formed links are:'
1310 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1311 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1317 if (ns.gt.0.and.dyn_ss) then
1321 forcon(i-nss)=forcon(i)
1328 dyn_ss_mask(iss(i))=.true.
1333 c write (iout,*) "DC"
1335 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1336 c & (dc(j,i+nres),j=1,3)
1338 if (i2ndstr.gt.0) call secstrp2dihc
1339 c call geom_to_var(nvar,x)
1340 c call etotal(energia(0))
1341 c call enerprint(energia(0))
1342 c call briefout(0,etot)
1344 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1345 cd write (iout,'(a)') 'Variable list:'
1346 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1348 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1349 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1350 & 'Processor',myrank,': end reading molecular data.'
1355 c--------------------------------------------------------------------------
1356 logical function seq_comp(itypea,itypeb,length)
1358 integer length,itypea(length),itypeb(length)
1361 if (itypea(i).ne.itypeb(i)) then
1369 c-----------------------------------------------------------------------------
1370 subroutine read_bridge
1371 C Read information about disulfide bridges.
1373 include 'DIMENSIONS'
1378 include 'COMMON.IOUNITS'
1379 include 'COMMON.GEO'
1380 include 'COMMON.VAR'
1381 include 'COMMON.INTERACT'
1382 include 'COMMON.LOCAL'
1383 include 'COMMON.NAMES'
1384 include 'COMMON.CHAIN'
1385 include 'COMMON.FFIELD'
1386 include 'COMMON.SBRIDGE'
1387 include 'COMMON.HEADER'
1388 include 'COMMON.CONTROL'
1389 include 'COMMON.DBASE'
1390 include 'COMMON.THREAD'
1391 include 'COMMON.TIME1'
1392 include 'COMMON.SETUP'
1394 C Read bridging residues.
1395 read (inp,*) ns,(iss(i),i=1,ns)
1397 if(me.eq.king.or..not.out1file)
1398 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1399 C Check whether the specified bridging residues are cystines.
1401 if (itype(iss(i)).ne.1) then
1402 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1403 & 'Do you REALLY think that the residue ',
1404 & restyp(itype(iss(i))),i,
1405 & ' can form a disulfide bridge?!!!'
1406 write (*,'(2a,i3,a)')
1407 & 'Do you REALLY think that the residue ',
1408 & restyp(itype(iss(i))),i,
1409 & ' can form a disulfide bridge?!!!'
1411 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1416 C Read preformed bridges.
1418 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1420 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1423 C Check if the residues involved in bridges are in the specified list of
1424 C bridging residues.
1427 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1428 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1429 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1430 & ' contains residues present in other pairs.'
1431 write (*,'(a,i3,a)') 'Disulfide pair',i,
1432 & ' contains residues present in other pairs.'
1434 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1440 if (ihpb(i).eq.iss(j)) goto 10
1442 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1445 if (jhpb(i).eq.iss(j)) goto 20
1447 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1453 ihpb(i)=ihpb(i)+nres
1454 jhpb(i)=jhpb(i)+nres
1460 c----------------------------------------------------------------------------
1461 subroutine read_x(kanal,*)
1463 include 'DIMENSIONS'
1464 include 'COMMON.GEO'
1465 include 'COMMON.VAR'
1466 include 'COMMON.CHAIN'
1467 include 'COMMON.IOUNITS'
1468 include 'COMMON.CONTROL'
1469 include 'COMMON.LOCAL'
1470 include 'COMMON.INTERACT'
1471 integer i,j,k,l,kanal
1472 c Read coordinates from input
1474 read(kanal,'(8f10.5)',end=10,err=10)
1475 & ((c(l,k),l=1,3),k=1,nres),
1476 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1479 c(j,2*nres)=c(j,nres)
1481 call int_from_cart1(.false.)
1484 dc(j,i)=c(j,i+1)-c(j,i)
1485 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1489 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1491 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1492 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1500 c----------------------------------------------------------------------------
1501 subroutine read_threadbase
1503 include 'DIMENSIONS'
1504 include 'COMMON.IOUNITS'
1505 include 'COMMON.GEO'
1506 include 'COMMON.VAR'
1507 include 'COMMON.INTERACT'
1508 include 'COMMON.LOCAL'
1509 include 'COMMON.NAMES'
1510 include 'COMMON.CHAIN'
1511 include 'COMMON.FFIELD'
1512 include 'COMMON.SBRIDGE'
1513 include 'COMMON.HEADER'
1514 include 'COMMON.CONTROL'
1515 include 'COMMON.DBASE'
1516 include 'COMMON.THREAD'
1517 include 'COMMON.TIME1'
1519 double precision dist
1520 C Read pattern database for threading.
1521 read (icbase,*) nseq
1523 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1524 & nres_base(2,i),nres_base(3,i)
1525 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1527 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1528 c & nres_base(2,i),nres_base(3,i)
1529 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1533 if (weidis.eq.0.0D0) weidis=0.1D0
1542 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1543 write (iout,'(a,i5)') 'nexcl: ',nexcl
1544 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1547 c------------------------------------------------------------------------------
1548 subroutine setup_var
1550 include 'DIMENSIONS'
1551 include 'COMMON.IOUNITS'
1552 include 'COMMON.GEO'
1553 include 'COMMON.VAR'
1554 include 'COMMON.INTERACT'
1555 include 'COMMON.LOCAL'
1556 include 'COMMON.NAMES'
1557 include 'COMMON.CHAIN'
1558 include 'COMMON.FFIELD'
1559 include 'COMMON.SBRIDGE'
1560 include 'COMMON.HEADER'
1561 include 'COMMON.CONTROL'
1562 include 'COMMON.DBASE'
1563 include 'COMMON.THREAD'
1564 include 'COMMON.TIME1'
1566 C Set up variable list.
1572 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1574 ialph(i,1)=nvar+nside
1578 if (indphi.gt.0) then
1580 else if (indback.gt.0) then
1587 c----------------------------------------------------------------------------
1588 subroutine gen_dist_constr
1589 C Generate CA distance constraints.
1591 include 'DIMENSIONS'
1592 include 'COMMON.IOUNITS'
1593 include 'COMMON.GEO'
1594 include 'COMMON.VAR'
1595 include 'COMMON.INTERACT'
1596 include 'COMMON.LOCAL'
1597 include 'COMMON.NAMES'
1598 include 'COMMON.CHAIN'
1599 include 'COMMON.FFIELD'
1600 include 'COMMON.SBRIDGE'
1601 include 'COMMON.HEADER'
1602 include 'COMMON.CONTROL'
1603 include 'COMMON.DBASE'
1604 include 'COMMON.THREAD'
1605 include 'COMMON.TIME1'
1606 integer i,j,itype_pdb(maxres)
1607 common /pizda/ itype_pdb
1608 double precision dist
1610 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1611 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1612 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1614 do i=nstart_sup,nstart_sup+nsup-1
1615 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1616 cd & ' seq_pdb', restyp(itype_pdb(i))
1617 do j=i+2,nstart_sup+nsup-1
1619 ihpb(nhpb)=i+nstart_seq-nstart_sup
1620 jhpb(nhpb)=j+nstart_seq-nstart_sup
1622 dhpb(nhpb)=dist(i,j)
1625 cd write (iout,'(a)') 'Distance constraints:'
1630 cd if (ii.gt.nres) then
1635 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1636 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1637 cd & dhpb(i),forcon(i)
1641 c----------------------------------------------------------------------------
1644 include 'DIMENSIONS'
1645 include 'COMMON.MAP'
1646 include 'COMMON.IOUNITS'
1648 character*3 angid(4) /'THE','PHI','ALP','OME'/
1649 character*80 mapcard,ucase
1651 read (inp,'(a)') mapcard
1652 mapcard=ucase(mapcard)
1653 if (index(mapcard,'PHI').gt.0) then
1655 else if (index(mapcard,'THE').gt.0) then
1657 else if (index(mapcard,'ALP').gt.0) then
1659 else if (index(mapcard,'OME').gt.0) then
1662 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1663 stop 'Error - illegal variable spec in MAP card.'
1665 call readi (mapcard,'RES1',res1(imap),0)
1666 call readi (mapcard,'RES2',res2(imap),0)
1667 if (res1(imap).eq.0) then
1668 res1(imap)=res2(imap)
1669 else if (res2(imap).eq.0) then
1670 res2(imap)=res1(imap)
1672 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1674 & 'Error - illegal definition of variable group in MAP.'
1675 stop 'Error - illegal definition of variable group in MAP.'
1677 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1678 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1679 call readi(mapcard,'NSTEP',nstep(imap),0)
1680 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1682 & 'Illegal boundary and/or step size specification in MAP.'
1683 stop 'Illegal boundary and/or step size specification in MAP.'
1688 c----------------------------------------------------------------------------
1691 include 'DIMENSIONS'
1692 include 'COMMON.IOUNITS'
1693 include 'COMMON.GEO'
1694 include 'COMMON.CSA'
1695 include 'COMMON.BANK'
1696 include 'COMMON.CONTROL'
1698 character*620 mcmcard
1699 call card_concat(mcmcard)
1701 call readi(mcmcard,'NCONF',nconf,50)
1702 call readi(mcmcard,'NADD',nadd,0)
1703 call readi(mcmcard,'JSTART',jstart,1)
1704 call readi(mcmcard,'JEND',jend,1)
1705 call readi(mcmcard,'NSTMAX',nstmax,500000)
1706 call readi(mcmcard,'N0',n0,1)
1707 call readi(mcmcard,'N1',n1,6)
1708 call readi(mcmcard,'N2',n2,4)
1709 call readi(mcmcard,'N3',n3,0)
1710 call readi(mcmcard,'N4',n4,0)
1711 call readi(mcmcard,'N5',n5,0)
1712 call readi(mcmcard,'N6',n6,10)
1713 call readi(mcmcard,'N7',n7,0)
1714 call readi(mcmcard,'N8',n8,0)
1715 call readi(mcmcard,'N9',n9,0)
1716 call readi(mcmcard,'N14',n14,0)
1717 call readi(mcmcard,'N15',n15,0)
1718 call readi(mcmcard,'N16',n16,0)
1719 call readi(mcmcard,'N17',n17,0)
1720 call readi(mcmcard,'N18',n18,0)
1722 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1724 call readi(mcmcard,'NDIFF',ndiff,2)
1725 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1726 call readi(mcmcard,'IS1',is1,1)
1727 call readi(mcmcard,'IS2',is2,8)
1728 call readi(mcmcard,'NRAN0',nran0,4)
1729 call readi(mcmcard,'NRAN1',nran1,2)
1730 call readi(mcmcard,'IRR',irr,1)
1731 call readi(mcmcard,'NSEED',nseed,20)
1732 call readi(mcmcard,'NTOTAL',ntotal,10000)
1733 call reada(mcmcard,'CUT1',cut1,2.0d0)
1734 call reada(mcmcard,'CUT2',cut2,5.0d0)
1735 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1736 call readi(mcmcard,'ICMAX',icmax,3)
1737 call readi(mcmcard,'IRESTART',irestart,0)
1738 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1741 call reada(mcmcard,'DELE',dele,20.0d0)
1742 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1743 call readi(mcmcard,'IREF',iref,0)
1744 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1745 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1746 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1747 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1748 write (iout,*) "NCONF_IN",nconf_in
1751 c----------------------------------------------------------------------------
1754 include 'DIMENSIONS'
1755 include 'COMMON.MCM'
1756 include 'COMMON.MCE'
1757 include 'COMMON.IOUNITS'
1759 character*320 mcmcard
1761 call card_concat(mcmcard)
1762 call readi(mcmcard,'MAXACC',maxacc,100)
1763 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1764 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1765 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1766 call readi(mcmcard,'MAXREPM',maxrepm,200)
1767 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1768 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1769 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1770 call reada(mcmcard,'E_UP',e_up,5.0D0)
1771 call reada(mcmcard,'DELTE',delte,0.1D0)
1772 call readi(mcmcard,'NSWEEP',nsweep,5)
1773 call readi(mcmcard,'NSTEPH',nsteph,0)
1774 call readi(mcmcard,'NSTEPC',nstepc,0)
1775 call reada(mcmcard,'TMIN',tmin,298.0D0)
1776 call reada(mcmcard,'TMAX',tmax,298.0D0)
1777 call readi(mcmcard,'NWINDOW',nwindow,0)
1778 call readi(mcmcard,'PRINT_MC',print_mc,0)
1779 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1780 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1781 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1782 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1783 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1784 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1785 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1786 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1787 if (nwindow.gt.0) then
1788 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1790 winlen(i)=winend(i)-winstart(i)+1
1793 if (tmax.lt.tmin) tmax=tmin
1794 if (tmax.eq.tmin) then
1798 if (nstepc.gt.0 .and. nsteph.gt.0) then
1799 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1800 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1802 C Probabilities of different move types
1803 sumpro_type(0)=0.0D0
1804 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1805 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1806 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1807 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1808 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1809 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1810 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1812 print *,'i',i,' sumprotype',sumpro_type(i)
1813 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1814 print *,'i',i,' sumprotype',sumpro_type(i)
1818 c----------------------------------------------------------------------------
1819 subroutine read_minim
1821 include 'DIMENSIONS'
1822 include 'COMMON.MINIM'
1823 include 'COMMON.IOUNITS'
1825 character*320 minimcard
1826 call card_concat(minimcard)
1827 call readi(minimcard,'MAXMIN',maxmin,2000)
1828 call readi(minimcard,'MAXFUN',maxfun,5000)
1829 call readi(minimcard,'MINMIN',minmin,maxmin)
1830 call readi(minimcard,'MINFUN',minfun,maxmin)
1831 call reada(minimcard,'TOLF',tolf,1.0D-2)
1832 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1833 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1834 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1835 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1836 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1837 & 'Options in energy minimization:'
1838 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1839 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1840 & 'MinMin:',MinMin,' MinFun:',MinFun,
1841 & ' TolF:',TolF,' RTolF:',RTolF
1844 c----------------------------------------------------------------------------
1845 subroutine read_angles(kanal,*)
1847 include 'DIMENSIONS'
1848 include 'COMMON.GEO'
1849 include 'COMMON.VAR'
1850 include 'COMMON.CHAIN'
1851 include 'COMMON.IOUNITS'
1852 include 'COMMON.CONTROL'
1854 c Read angles from input
1856 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1857 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1858 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1859 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1862 c 9/7/01 avoid 180 deg valence angle
1863 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1865 theta(i)=deg2rad*theta(i)
1866 phi(i)=deg2rad*phi(i)
1867 alph(i)=deg2rad*alph(i)
1868 omeg(i)=deg2rad*omeg(i)
1873 c----------------------------------------------------------------------------
1874 subroutine reada(rekord,lancuch,wartosc,default)
1876 character*(*) rekord,lancuch
1877 double precision wartosc,default
1880 iread=index(rekord,lancuch)
1881 if (iread.eq.0) then
1885 iread=iread+ilen(lancuch)+1
1886 read (rekord(iread:),*,err=10,end=10) wartosc
1891 c----------------------------------------------------------------------------
1892 subroutine readi(rekord,lancuch,wartosc,default)
1894 character*(*) rekord,lancuch
1895 integer wartosc,default
1898 iread=index(rekord,lancuch)
1899 if (iread.eq.0) then
1903 iread=iread+ilen(lancuch)+1
1904 read (rekord(iread:),*,err=10,end=10) wartosc
1909 c----------------------------------------------------------------------------
1910 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1913 integer tablica(dim),default
1914 character*(*) rekord,lancuch
1921 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1922 if (iread.eq.0) return
1923 iread=iread+ilen(lancuch)+1
1924 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1927 c----------------------------------------------------------------------------
1928 subroutine multreada(rekord,lancuch,tablica,dim,default)
1931 double precision tablica(dim),default
1932 character*(*) rekord,lancuch
1939 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1940 if (iread.eq.0) return
1941 iread=iread+ilen(lancuch)+1
1942 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1945 c----------------------------------------------------------------------------
1946 subroutine openunits
1948 include 'DIMENSIONS'
1952 character*16 form,nodename
1955 include 'COMMON.SETUP'
1956 include 'COMMON.IOUNITS'
1958 include 'COMMON.CONTROL'
1959 integer lenpre,lenpot,ilen,lentmp,npos
1961 character*3 out1file_text,ucase
1966 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1967 call getenv_loc("PREFIX",prefix)
1969 call getenv_loc("POT",pot)
1970 call getenv_loc("DIRTMP",tmpdir)
1971 call getenv_loc("CURDIR",curdir)
1972 call getenv_loc("OUT1FILE",out1file_text)
1973 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1974 out1file_text=ucase(out1file_text)
1975 if (out1file_text(1:1).eq."Y") then
1978 out1file=fg_rank.gt.0
1983 if (lentmp.gt.0) then
1984 write (*,'(80(1h!))')
1985 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1986 write (*,'(80(1h!))')
1987 write (*,*)"All output files will be on node /tmp directory."
1989 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1990 if (me.eq.king) then
1991 write (*,*) "The master node is ",nodename
1992 else if (fg_rank.eq.0) then
1993 write (*,*) "I am the CG slave node ",nodename
1995 write (*,*) "I am the FG slave node ",nodename
1998 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1999 lenpre = lentmp+lenpre+1
2001 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2002 C Get the names and open the input files
2003 #if defined(WINIFL) || defined(WINPGI)
2004 open(1,file=pref_orig(:ilen(pref_orig))//
2005 & '.inp',status='old',readonly,shared)
2006 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2007 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2008 C Get parameter filenames and open the parameter files.
2009 call getenv_loc('BONDPAR',bondname)
2010 open (ibond,file=bondname,status='old',readonly,shared)
2011 call getenv_loc('THETPAR',thetname)
2012 open (ithep,file=thetname,status='old',readonly,shared)
2013 call getenv_loc('ROTPAR',rotname)
2014 open (irotam,file=rotname,status='old',readonly,shared)
2015 call getenv_loc('TORPAR',torname)
2016 open (itorp,file=torname,status='old',readonly,shared)
2017 call getenv_loc('TORDPAR',tordname)
2018 open (itordp,file=tordname,status='old',readonly,shared)
2019 call getenv_loc('FOURIER',fouriername)
2020 open (ifourier,file=fouriername,status='old',readonly,shared)
2021 call getenv_loc('ELEPAR',elename)
2022 open (ielep,file=elename,status='old',readonly,shared)
2023 call getenv_loc('SIDEPAR',sidename)
2024 open (isidep,file=sidename,status='old',readonly,shared)
2025 call getenv_loc('LIPTRANPAR',liptranname)
2026 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2027 #elif (defined CRAY) || (defined AIX)
2028 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2030 c print *,"Processor",myrank," opened file 1"
2031 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2032 c print *,"Processor",myrank," opened file 9"
2033 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2034 C Get parameter filenames and open the parameter files.
2035 call getenv_loc('BONDPAR',bondname)
2036 open (ibond,file=bondname,status='old',action='read')
2037 c print *,"Processor",myrank," opened file IBOND"
2038 call getenv_loc('THETPAR',thetname)
2039 open (ithep,file=thetname,status='old',action='read')
2041 call getenv_loc('THETPARPDB',thetname_pdb)
2042 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2044 c print *,"Processor",myrank," opened file ITHEP"
2045 call getenv_loc('ROTPAR',rotname)
2046 open (irotam,file=rotname,status='old',action='read')
2048 call getenv_loc('ROTPARPDB',rotname_pdb)
2049 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2051 c print *,"Processor",myrank," opened file IROTAM"
2052 call getenv_loc('TORPAR',torname)
2053 open (itorp,file=torname,status='old',action='read')
2054 c print *,"Processor",myrank," opened file ITORP"
2055 call getenv_loc('TORDPAR',tordname)
2056 open (itordp,file=tordname,status='old',action='read')
2057 c print *,"Processor",myrank," opened file ITORDP"
2058 call getenv_loc('SCCORPAR',sccorname)
2059 open (isccor,file=sccorname,status='old',action='read')
2060 c print *,"Processor",myrank," opened file ISCCOR"
2061 call getenv_loc('FOURIER',fouriername)
2062 open (ifourier,file=fouriername,status='old',action='read')
2063 c print *,"Processor",myrank," opened file IFOURIER"
2064 call getenv_loc('ELEPAR',elename)
2065 open (ielep,file=elename,status='old',action='read')
2066 c print *,"Processor",myrank," opened file IELEP"
2067 call getenv_loc('SIDEPAR',sidename)
2068 open (isidep,file=sidename,status='old',action='read')
2069 call getenv_loc('LIPTRANPAR',liptranname)
2070 open (iliptranpar,file=liptranname,status='old',action='read')
2071 c print *,"Processor",myrank," opened file ISIDEP"
2072 c print *,"Processor",myrank," opened parameter files"
2074 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2075 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2076 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2077 C Get parameter filenames and open the parameter files.
2078 call getenv_loc('BONDPAR',bondname)
2079 open (ibond,file=bondname,status='old')
2080 call getenv_loc('THETPAR',thetname)
2081 open (ithep,file=thetname,status='old')
2083 call getenv_loc('THETPARPDB',thetname_pdb)
2084 open (ithep_pdb,file=thetname_pdb,status='old')
2086 call getenv_loc('ROTPAR',rotname)
2087 open (irotam,file=rotname,status='old')
2089 call getenv_loc('ROTPARPDB',rotname_pdb)
2090 open (irotam_pdb,file=rotname_pdb,status='old')
2092 call getenv_loc('TORPAR',torname)
2093 open (itorp,file=torname,status='old')
2094 call getenv_loc('TORDPAR',tordname)
2095 open (itordp,file=tordname,status='old')
2096 call getenv_loc('SCCORPAR',sccorname)
2097 open (isccor,file=sccorname,status='old')
2098 call getenv_loc('FOURIER',fouriername)
2099 open (ifourier,file=fouriername,status='old')
2100 call getenv_loc('ELEPAR',elename)
2101 open (ielep,file=elename,status='old')
2102 call getenv_loc('SIDEPAR',sidename)
2103 open (isidep,file=sidename,status='old')
2104 call getenv_loc('LIPTRANPAR',liptranname)
2105 open (iliptranpar,file=liptranname,status='old')
2107 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2109 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
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',readonly)
2114 call getenv_loc('THETPAR',thetname)
2115 open (ithep,file=thetname,status='old',readonly)
2116 call getenv_loc('ROTPAR',rotname)
2117 open (irotam,file=rotname,status='old',readonly)
2118 call getenv_loc('TORPAR',torname)
2119 open (itorp,file=torname,status='old',readonly)
2120 call getenv_loc('TORDPAR',tordname)
2121 open (itordp,file=tordname,status='old',readonly)
2122 call getenv_loc('SCCORPAR',sccorname)
2123 open (isccor,file=sccorname,status='old',readonly)
2125 call getenv_loc('THETPARPDB',thetname_pdb)
2126 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2128 call getenv_loc('FOURIER',fouriername)
2129 open (ifourier,file=fouriername,status='old',readonly)
2130 call getenv_loc('ELEPAR',elename)
2131 open (ielep,file=elename,status='old',readonly)
2132 call getenv_loc('SIDEPAR',sidename)
2133 open (isidep,file=sidename,status='old',readonly)
2134 call getenv_loc('LIPTRANPAR',liptranname)
2135 open (iliptranpar,file=liptranname,status='old',action='read')
2137 call getenv_loc('ROTPARPDB',rotname_pdb)
2138 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2142 call getenv_loc('TUBEPAR',tubename)
2143 #if defined(WINIFL) || defined(WINPGI)
2144 open (itube,file=tubename,status='old',readonly,shared)
2145 #elif (defined CRAY) || (defined AIX)
2146 open (itube,file=tubename,status='old',action='read')
2148 open (itube,file=tubename,status='old')
2150 open (itube,file=tubename,status='old',readonly)
2155 C 8/9/01 In the newest version SCp interaction constants are read from a file
2156 C Use -DOLDSCP to use hard-coded constants instead.
2158 call getenv_loc('SCPPAR',scpname)
2159 #if defined(WINIFL) || defined(WINPGI)
2160 open (iscpp,file=scpname,status='old',readonly,shared)
2161 #elif (defined CRAY) || (defined AIX)
2162 open (iscpp,file=scpname,status='old',action='read')
2164 open (iscpp,file=scpname,status='old')
2166 open (iscpp,file=scpname,status='old',readonly)
2169 call getenv_loc('PATTERN',patname)
2170 #if defined(WINIFL) || defined(WINPGI)
2171 open (icbase,file=patname,status='old',readonly,shared)
2172 #elif (defined CRAY) || (defined AIX)
2173 open (icbase,file=patname,status='old',action='read')
2175 open (icbase,file=patname,status='old')
2177 open (icbase,file=patname,status='old',readonly)
2180 C Open output file only for CG processes
2181 c print *,"Processor",myrank," fg_rank",fg_rank
2182 if (fg_rank.eq.0) then
2184 if (nodes.eq.1) then
2187 npos = dlog10(dfloat(nodes-1))+1
2189 if (npos.lt.3) npos=3
2190 write (liczba,'(i1)') npos
2191 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2193 write (liczba,form) me
2194 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2195 & liczba(:ilen(liczba))
2196 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2198 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2200 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2201 & liczba(:ilen(liczba))//'.mol2'
2202 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2203 & liczba(:ilen(liczba))//'.stat'
2205 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2206 & //liczba(:ilen(liczba))//'.stat')
2207 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2210 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2211 & liczba(:ilen(liczba))//'.const'
2216 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2217 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2218 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2219 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2220 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2222 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2224 rest2name=prefix(:ilen(prefix))//'.rst'
2226 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2229 #if defined(AIX) || defined(PGI) || defined(CRAY)
2230 if (me.eq.king .or. .not. out1file) then
2231 open(iout,file=outname,status='unknown')
2233 open(iout,file="/dev/null",status="unknown")
2237 if (fg_rank.gt.0) then
2238 write (liczba,'(i3.3)') myrank/nfgtasks
2239 write (ll,'(bz,i3.3)') fg_rank
2240 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2246 open(igeom,file=intname,status='unknown',position='append')
2247 open(ipdb,file=pdbname,status='unknown')
2248 open(imol2,file=mol2name,status='unknown')
2249 open(istat,file=statname,status='unknown',position='append')
2251 c1out open(iout,file=outname,status='unknown')
2254 if (me.eq.king .or. .not.out1file)
2255 & open(iout,file=outname,status='unknown')
2257 if (fg_rank.gt.0) then
2258 write (liczba,'(i3.3)') myrank/nfgtasks
2259 write (ll,'(bz,i3.3)') fg_rank
2260 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2265 open(igeom,file=intname,status='unknown',access='append')
2266 open(ipdb,file=pdbname,status='unknown')
2267 open(imol2,file=mol2name,status='unknown')
2268 open(istat,file=statname,status='unknown',access='append')
2270 c1out open(iout,file=outname,status='unknown')
2273 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2274 csa_seed=prefix(:lenpre)//'.CSA.seed'
2275 csa_history=prefix(:lenpre)//'.CSA.history'
2276 csa_bank=prefix(:lenpre)//'.CSA.bank'
2277 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2278 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2279 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2280 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2281 csa_int=prefix(:lenpre)//'.int'
2282 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2283 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2284 csa_in=prefix(:lenpre)//'.CSA.in'
2285 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2288 write (iout,'(80(1h-))')
2289 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2290 write (iout,'(80(1h-))')
2291 write (iout,*) "Input file : ",
2292 & pref_orig(:ilen(pref_orig))//'.inp'
2293 write (iout,*) "Output file : ",
2294 & outname(:ilen(outname))
2296 write (iout,*) "Sidechain potential file : ",
2297 & sidename(:ilen(sidename))
2299 write (iout,*) "SCp potential file : ",
2300 & scpname(:ilen(scpname))
2302 write (iout,*) "Electrostatic potential file : ",
2303 & elename(:ilen(elename))
2304 write (iout,*) "Cumulant coefficient file : ",
2305 & fouriername(:ilen(fouriername))
2306 write (iout,*) "Torsional parameter file : ",
2307 & torname(:ilen(torname))
2308 write (iout,*) "Double torsional parameter file : ",
2309 & tordname(:ilen(tordname))
2310 write (iout,*) "SCCOR parameter file : ",
2311 & sccorname(:ilen(sccorname))
2312 write (iout,*) "Bond & inertia constant file : ",
2313 & bondname(:ilen(bondname))
2314 write (iout,*) "Bending-torsion parameter file : ",
2315 & thetname(:ilen(thetname))
2316 write (iout,*) "Rotamer parameter file : ",
2317 & rotname(:ilen(rotname))
2318 write (iout,*) "Thetpdb parameter file : ",
2319 & thetname_pdb(:ilen(thetname_pdb))
2320 write (iout,*) "Threading database : ",
2321 & patname(:ilen(patname))
2323 &write (iout,*)" DIRTMP : ",
2325 write (iout,'(80(1h-))')
2329 c----------------------------------------------------------------------------
2330 subroutine card_concat(card)
2331 implicit real*8 (a-h,o-z)
2332 include 'DIMENSIONS'
2333 include 'COMMON.IOUNITS'
2335 character*80 karta,ucase
2337 read (inp,'(a)') karta
2340 do while (karta(80:80).eq.'&')
2341 card=card(:ilen(card)+1)//karta(:79)
2342 read (inp,'(a)') karta
2345 card=card(:ilen(card)+1)//karta
2348 c------------------------------------------------------------------------------
2351 include 'DIMENSIONS'
2352 include 'COMMON.CHAIN'
2353 include 'COMMON.IOUNITS'
2354 include 'COMMON.CONTROL'
2356 include 'COMMON.QRESTR'
2358 include 'COMMON.LAGRANGE.5diag'
2360 include 'COMMON.LAGRANGE'
2363 open(irest2,file=rest2name,status='unknown')
2364 read(irest2,*) totT,EK,potE,totE,t_bath
2367 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2370 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2373 read (irest2,*) iset
2378 c------------------------------------------------------------------------------
2379 subroutine read_fragments
2381 include 'DIMENSIONS'
2385 include 'COMMON.SETUP'
2386 include 'COMMON.CHAIN'
2387 include 'COMMON.IOUNITS'
2389 include 'COMMON.QRESTR'
2390 include 'COMMON.CONTROL'
2392 read(inp,*) nset,nfrag,npair,nfrag_back
2393 loc_qlike=(nfrag_back.lt.0)
2394 nfrag_back=iabs(nfrag_back)
2395 if(me.eq.king.or..not.out1file)
2396 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2397 & " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2399 read(inp,*) mset(iset)
2401 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2403 if(me.eq.king.or..not.out1file)
2404 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2405 & ifrag(2,i,iset), qinfrag(i,iset)
2408 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2410 if(me.eq.king.or..not.out1file)
2411 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2412 & ipair(2,i,iset), qinpair(i,iset)
2416 read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2417 & wfrag_back(2,i,iset),qin_back(2,i,iset),
2418 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2419 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2420 if(me.eq.king.or..not.out1file)
2421 & write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2422 & wfrag_back(2,i,iset),qin_back(3,i,iset),
2423 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2424 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2428 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2429 & wfrag_back(3,i,iset),
2430 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2431 if(me.eq.king.or..not.out1file)
2432 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2433 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2439 C---------------------------------------------------------------------------
2440 subroutine read_afminp
2442 include 'DIMENSIONS'
2446 include 'COMMON.SETUP'
2447 include 'COMMON.CONTROL'
2448 include 'COMMON.CHAIN'
2449 include 'COMMON.IOUNITS'
2450 include 'COMMON.SBRIDGE'
2451 character*320 afmcard
2453 c print *, "wchodze"
2454 call card_concat(afmcard)
2455 call readi(afmcard,"BEG",afmbeg,0)
2456 call readi(afmcard,"END",afmend,0)
2457 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2458 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2459 print *,'FORCE=' ,forceAFMconst
2460 CCCC NOW PROPERTIES FOR AFM
2463 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2465 distafminit=dsqrt(distafminit)
2466 c print *,'initdist',distafminit
2469 c-------------------------------------------------------------------------------
2470 subroutine read_saxs_constr
2472 include 'DIMENSIONS'
2476 include 'COMMON.SETUP'
2477 include 'COMMON.CONTROL'
2478 include 'COMMON.SAXS'
2479 include 'COMMON.CHAIN'
2480 include 'COMMON.IOUNITS'
2481 include 'COMMON.SBRIDGE'
2482 double precision cm(3),cnorm
2485 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2487 if (saxs_mode.eq.0) then
2488 c SAXS distance distribution
2490 read(inp,*) distsaxs(i),Psaxs(i)
2494 Cnorm = Cnorm + Psaxs(i)
2496 write (iout,*) "Cnorm",Cnorm
2498 Psaxs(i)=Psaxs(i)/Cnorm
2500 write (iout,*) "Normalized distance distribution from SAXS"
2502 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2506 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2508 write (iout,*) "Wsaxs0",Wsaxs0
2512 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2519 cm(j)=cm(j)+Csaxs(j,i)
2527 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2530 write (iout,*) "SAXS sphere coordinates"
2532 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2538 c-------------------------------------------------------------------------------
2539 subroutine read_dist_constr
2541 include 'DIMENSIONS'
2545 include 'COMMON.SETUP'
2546 include 'COMMON.CONTROL'
2547 include 'COMMON.CHAIN'
2548 include 'COMMON.IOUNITS'
2549 include 'COMMON.SBRIDGE'
2550 include 'COMMON.INTERACT'
2551 integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
2552 integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
2553 double precision wfrag_(100),wpair_(1000)
2554 double precision ddjk,dist,dist_cut,fordepthmax
2555 character*5000 controlcard
2556 logical normalize,next
2558 double precision scal_bfac
2559 double precision xlink(4,0:4) /
2561 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2562 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2563 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2564 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2565 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2566 c print *, "WCHODZE"
2567 write (iout,*) "Calling read_dist_constr"
2568 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2570 restr_on_coord=.false.
2579 call card_concat(controlcard)
2580 next = index(controlcard,"NEXT").gt.0
2581 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2582 write (iout,*) "restr_type",restr_type
2583 call readi(controlcard,"NFRAG",nfrag_,0)
2584 call readi(controlcard,"NFRAG",nfrag_,0)
2585 call readi(controlcard,"NPAIR",npair_,0)
2586 call readi(controlcard,"NDIST",ndist_,0)
2587 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2588 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2589 if (restr_type.eq.10)
2590 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2591 if (restr_type.eq.12)
2592 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2593 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2594 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2595 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2596 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2597 normalize = index(controlcard,"NORMALIZE").gt.0
2598 write (iout,*) "WBOLTZD",wboltzd
2599 write (iout,*) "SCAL_PEAK",scal_peak
2600 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2601 write (iout,*) "IFRAG"
2603 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2605 write (iout,*) "IPAIR"
2607 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2609 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2611 & "Distance restraints as generated from reference structure"
2613 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2614 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2615 & ifrag_(2,i)=nstart_sup+nsup-1
2616 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2618 if (wfrag_(i).eq.0.0d0) cycle
2619 do j=ifrag_(1,i),ifrag_(2,i)-1
2620 do k=j+1,ifrag_(2,i)
2621 c write (iout,*) "j",j," k",k
2623 if (restr_type.eq.1) then
2629 forcon(nhpb)=wfrag_(i)
2630 else if (constr_dist.eq.2) then
2631 if (ddjk.le.dist_cut) then
2637 forcon(nhpb)=wfrag_(i)
2639 else if (restr_type.eq.3) then
2645 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2648 if (.not.out1file .or. me.eq.king)
2649 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2650 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2652 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2653 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2659 if (wpair_(i).eq.0.0d0) cycle
2667 do j=ifrag_(1,ii),ifrag_(2,ii)
2668 do k=ifrag_(1,jj),ifrag_(2,jj)
2670 if (restr_type.eq.1) then
2676 forcon(nhpb)=wpair_(i)
2677 else if (constr_dist.eq.2) then
2678 if (ddjk.le.dist_cut) then
2684 forcon(nhpb)=wpair_(i)
2686 else if (restr_type.eq.3) then
2692 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2695 if (.not.out1file .or. me.eq.king)
2696 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2697 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2699 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2700 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2707 write (iout,*) "Distance restraints as read from input"
2709 if (restr_type.eq.12) then
2710 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2711 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2712 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2713 & fordepth_peak(nhpb_peak+1),npeak
2714 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2715 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2716 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2717 c & fordepth_peak(nhpb_peak+1),npeak
2718 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2719 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2720 nhpb_peak=nhpb_peak+1
2721 irestr_type_peak(nhpb_peak)=12
2722 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2725 if (.not.out1file .or. me.eq.king)
2726 & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2727 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2728 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2729 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2730 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2732 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2733 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2734 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2735 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2736 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2738 if (ibecarb_peak(nhpb_peak).eq.3) then
2739 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2740 else if (ibecarb_peak(nhpb_peak).eq.2) then
2741 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2742 else if (ibecarb_peak(nhpb_peak).eq.1) then
2743 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2744 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2746 else if (restr_type.eq.11) then
2747 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2748 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2749 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2750 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2752 irestr_type(nhpb)=11
2754 if (.not.out1file .or. me.eq.king)
2755 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2756 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2757 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2759 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2760 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2761 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2763 c if (ibecarb(nhpb).gt.0) then
2764 c ihpb(nhpb)=ihpb(nhpb)+nres
2765 c jhpb(nhpb)=jhpb(nhpb)+nres
2767 if (ibecarb(nhpb).eq.3) then
2768 ihpb(nhpb)=ihpb(nhpb)+nres
2769 else if (ibecarb(nhpb).eq.2) then
2770 ihpb(nhpb)=ihpb(nhpb)+nres
2771 else if (ibecarb(nhpb).eq.1) then
2772 ihpb(nhpb)=ihpb(nhpb)+nres
2773 jhpb(nhpb)=jhpb(nhpb)+nres
2775 else if (restr_type.eq.10) then
2776 c Cross-lonk Markov-like potential
2777 call card_concat(controlcard)
2778 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2779 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2781 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2782 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2783 if (index(controlcard,"ZL").gt.0) then
2785 else if (index(controlcard,"ADH").gt.0) then
2787 else if (index(controlcard,"PDH").gt.0) then
2789 else if (index(controlcard,"DSS").gt.0) then
2794 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2795 & xlink(1,link_type))
2796 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2797 & xlink(2,link_type))
2798 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2799 & xlink(3,link_type))
2800 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2801 & xlink(4,link_type))
2802 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2803 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2804 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2805 if (forcon(nhpb+1).le.0.0d0 .or.
2806 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2808 irestr_type(nhpb)=10
2809 if (ibecarb(nhpb).eq.3) then
2810 jhpb(nhpb)=jhpb(nhpb)+nres
2811 else if (ibecarb(nhpb).eq.2) then
2812 ihpb(nhpb)=ihpb(nhpb)+nres
2813 else if (ibecarb(nhpb).eq.1) then
2814 ihpb(nhpb)=ihpb(nhpb)+nres
2815 jhpb(nhpb)=jhpb(nhpb)+nres
2818 if (.not.out1file .or. me.eq.king)
2819 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2820 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2821 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2824 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2825 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2826 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2831 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2832 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2833 if (forcon(nhpb+1).gt.0.0d0) then
2835 if (dhpb1(nhpb).eq.0.0d0) then
2840 if (ibecarb(nhpb).eq.3) then
2841 jhpb(nhpb)=jhpb(nhpb)+nres
2842 else if (ibecarb(nhpb).eq.2) then
2843 ihpb(nhpb)=ihpb(nhpb)+nres
2844 else if (ibecarb(nhpb).eq.1) then
2845 ihpb(nhpb)=ihpb(nhpb)+nres
2846 jhpb(nhpb)=jhpb(nhpb)+nres
2848 if (dhpb(nhpb).eq.0.0d0)
2849 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2852 if (.not.out1file .or. me.eq.king)
2853 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2854 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2856 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2857 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2860 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2861 C if (forcon(nhpb+1).gt.0.0d0) then
2863 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2866 if (restr_type.eq.4) then
2867 write (iout,*) "The BFAC array"
2869 write (iout,'(i5,f10.5)') i,bfac(i)
2872 if (itype(i).eq.ntyp1) cycle
2874 if (itype(j).eq.ntyp1) cycle
2875 if (itype(i).eq.10) then
2880 if (itype(j).eq.10) then
2890 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2893 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2894 ihpb(nhpb)=i+nres*ii
2895 jhpb(nhpb)=j+nres*jj
2896 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2898 if (.not.out1file .or. me.eq.king) then
2899 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2900 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2901 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2905 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2906 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2907 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2917 if (restr_type.eq.5) then
2918 restr_on_coord=.true.
2920 if (itype(i).eq.ntyp1) cycle
2921 bfac(i)=(scal_bfac/bfac(i))**2
2930 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2931 & fordepthmax=fordepth(i)
2934 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
2939 c-------------------------------------------------------------------------------
2940 subroutine read_constr_homology
2942 include 'DIMENSIONS'
2946 include 'COMMON.SETUP'
2947 include 'COMMON.CONTROL'
2948 include 'COMMON.HOMOLOGY'
2949 include 'COMMON.CHAIN'
2950 include 'COMMON.IOUNITS'
2952 include 'COMMON.QRESTR'
2953 include 'COMMON.GEO'
2954 include 'COMMON.INTERACT'
2955 include 'COMMON.NAMES'
2957 c For new homol impl
2959 include 'COMMON.VAR'
2962 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2964 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2965 c & sigma_odl_temp(maxres,maxres,max_template)
2967 character*24 model_ki_dist, model_ki_angle
2968 character*500 controlcard
2969 integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
2970 & ik,iistart,iishift
2975 c FP - Nov. 2014 Temporary specifications for new vars
2977 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2979 double precision, dimension (max_template,maxres) :: rescore
2980 double precision, dimension (max_template,maxres) :: rescore2
2981 double precision, dimension (max_template,maxres) :: rescore3
2982 double precision distal
2983 character*24 pdbfile,tpl_k_rescore
2984 c -----------------------------------------------------------------
2985 c Reading multiple PDB ref structures and calculation of retraints
2986 c not using pre-computed ones stored in files model_ki_{dist,angle}
2988 c -----------------------------------------------------------------
2991 c Alternative: reading from input
2992 call card_concat(controlcard)
2993 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2994 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2995 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2996 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2997 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2998 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2999 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3000 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3001 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3002 if(.not.read2sigma.and.start_from_model) then
3003 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3004 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3005 start_from_model=.false.
3007 if(start_from_model .and. (me.eq.king .or. .not. out1file))
3008 & write(iout,*) 'START_FROM_MODELS is ON'
3009 if(start_from_model .and. rest) then
3010 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3011 write(iout,*) 'START_FROM_MODELS is OFF'
3012 write(iout,*) 'remove restart keyword from input'
3015 if (homol_nset.gt.1)then
3016 call card_concat(controlcard)
3017 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3018 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3019 write(iout,*) "iset homology_weight "
3021 write(iout,*) i,waga_homology(i)
3024 iset=mod(kolor,homol_nset)+1
3027 waga_homology(1)=1.0
3030 cd write (iout,*) "nnt",nnt," nct",nct
3037 c write(iout,*) 'nnt=',nnt,'nct=',nct
3040 do k=1,constr_homology
3053 if (read_homol_frag) then
3054 call read_klapaucjusz
3057 do k=1,constr_homology
3059 read(inp,'(a)') pdbfile
3060 if(me.eq.king .or. .not. out1file)
3061 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3062 & pdbfile(:ilen(pdbfile))
3063 open(ipdbin,file=pdbfile,status='old',err=33)
3065 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3066 & pdbfile(:ilen(pdbfile))
3069 c print *,'Begin reading pdb data'
3071 c Files containing res sim or local scores (former containing sigmas)
3074 write(kic2,'(bz,i2.2)') k
3076 tpl_k_rescore="template"//kic2//".sco"
3079 if (read2sigma) then
3080 call readpdb_template(k)
3085 c Distance restraints
3088 C Copy the coordinates from reference coordinates (?)
3092 c write (iout,*) "c(",j,i,") =",c(j,i)
3096 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3098 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3099 open (ientin,file=tpl_k_rescore,status='old')
3100 if (nnt.gt.1) rescore(k,1)=0.0d0
3101 do irec=nnt,nct ! loop for reading res sim
3102 if (read2sigma) then
3103 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3104 & rescore3_tmp,idomain_tmp
3106 idomain(k,i_tmp)=idomain_tmp
3107 rescore(k,i_tmp)=rescore_tmp
3108 rescore2(k,i_tmp)=rescore2_tmp
3109 rescore3(k,i_tmp)=rescore3_tmp
3110 if (.not. out1file .or. me.eq.king)
3111 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3112 & i_tmp,rescore2_tmp,rescore_tmp,
3113 & rescore3_tmp,idomain_tmp
3116 read (ientin,*,end=1401) rescore_tmp
3118 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3119 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3120 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3125 if (waga_dist.ne.0.0d0) then
3133 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3134 c write (iout,*) k,i,j,distal,dist2_cut
3136 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3137 & .and. distal.le.dist2_cut ) then
3143 c write (iout,*) "k",k
3144 c write (iout,*) "i",i," j",j," constr_homology",
3149 if (read2sigma) then
3152 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3154 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3155 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3156 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3158 if (odl(k,ii).le.dist_cut) then
3159 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3162 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3163 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3165 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3166 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3170 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3173 l_homo(k,ii)=.false.
3180 c Theta, dihedral and SC retraints
3182 if (waga_angle.gt.0.0d0) then
3183 c open (ientin,file=tpl_k_sigma_dih,status='old')
3184 c do irec=1,maxres-3 ! loop for reading sigma_dih
3185 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3186 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3187 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3188 c & sigma_dih(k,i+nnt-1)
3193 if (idomain(k,i).eq.0) then
3197 dih(k,i)=phiref(i) ! right?
3198 c read (ientin,*) sigma_dih(k,i) ! original variant
3199 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3200 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3201 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3202 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3203 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3205 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3206 & rescore(k,i-2)+rescore(k,i-3))/4.0
3207 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3208 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3209 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3210 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3211 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3212 c Instead of res sim other local measure of b/b str reliability possible
3213 if (sigma_dih(k,i).ne.0)
3214 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3215 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3220 if (waga_theta.gt.0.0d0) then
3221 c open (ientin,file=tpl_k_sigma_theta,status='old')
3222 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3223 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3224 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3225 c & sigma_theta(k,i+nnt-1)
3230 do i = nnt+2,nct ! right? without parallel.
3231 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3232 c do i=ithet_start,ithet_end ! with FG parallel.
3233 if (idomain(k,i).eq.0) then
3234 sigma_theta(k,i)=0.0
3237 thetatpl(k,i)=thetaref(i)
3238 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3239 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3240 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3241 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3242 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3243 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3244 & rescore(k,i-2))/3.0
3245 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3246 if (sigma_theta(k,i).ne.0)
3247 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3249 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3250 c rescore(k,i-2) ! right expression ?
3251 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3255 if (waga_d.gt.0.0d0) then
3256 c open (ientin,file=tpl_k_sigma_d,status='old')
3257 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3258 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3259 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3260 c & sigma_d(k,i+nnt-1)
3264 do i = nnt,nct ! right? without parallel.
3265 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3266 c do i=loc_start,loc_end ! with FG parallel.
3267 if (itype(i).eq.10) cycle
3268 if (idomain(k,i).eq.0 ) then
3275 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3276 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3277 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3278 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3279 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3280 if (sigma_d(k,i).ne.0)
3281 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3283 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3284 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3285 c read (ientin,*) sigma_d(k,i) ! 1st variant
3290 c remove distance restraints not used in any model from the list
3291 c shift data in all arrays
3293 if (waga_dist.ne.0.0d0) then
3299 if (ii_in_use(ii).eq.0.and.liiflag) then
3303 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3304 & .not.liiflag.and.ii.eq.lim_odl) then
3305 if (ii.eq.lim_odl) then
3306 iishift=ii-iistart+1
3311 do ki=iistart,lim_odl-iishift
3312 ires_homo(ki)=ires_homo(ki+iishift)
3313 jres_homo(ki)=jres_homo(ki+iishift)
3314 ii_in_use(ki)=ii_in_use(ki+iishift)
3315 do k=1,constr_homology
3316 odl(k,ki)=odl(k,ki+iishift)
3317 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3318 l_homo(k,ki)=l_homo(k,ki+iishift)
3322 lim_odl=lim_odl-iishift
3328 endif ! .not. klapaucjusz
3330 if (constr_homology.gt.0) call homology_partition
3331 if (constr_homology.gt.0) call init_int_table
3332 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3333 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3337 if (.not.out_template_restr) return
3338 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3339 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3340 write (iout,*) "Distance restraints from templates"
3342 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3343 & ii,ires_homo(ii),jres_homo(ii),
3344 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3345 & ki=1,constr_homology)
3347 write (iout,*) "Dihedral angle restraints from templates"
3349 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3350 & (rad2deg*dih(ki,i),
3351 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3353 write (iout,*) "Virtual-bond angle restraints from templates"
3355 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3356 & (rad2deg*thetatpl(ki,i),
3357 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3359 write (iout,*) "SC restraints from templates"
3361 write(iout,'(i5,100(4f8.2,4x))') i,
3362 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3363 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3366 c -----------------------------------------------------------------
3369 c----------------------------------------------------------------------
3371 subroutine flush(iu)
3376 subroutine flush(iu)
3381 c------------------------------------------------------------------------------
3382 subroutine copy_to_tmp(source)
3384 include "DIMENSIONS"
3385 include "COMMON.IOUNITS"
3386 character*(*) source
3387 character* 256 tmpfile
3391 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3392 inquire(file=tmpfile,exist=ex)
3394 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3395 & " to temporary directory..."
3396 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3397 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3401 c------------------------------------------------------------------------------
3402 subroutine move_from_tmp(source)
3404 include "DIMENSIONS"
3405 include "COMMON.IOUNITS"
3406 character*(*) source
3409 write (*,*) "Moving ",source(:ilen(source)),
3410 & " from temporary directory to working directory"
3411 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3412 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3415 c------------------------------------------------------------------------------
3416 subroutine random_init(seed)
3418 C Initialize random number generator
3421 include 'DIMENSIONS'
3424 logical OKRandom, prng_restart
3426 integer iseed_array(4)
3427 integer error_msg,ierr
3429 include 'COMMON.IOUNITS'
3430 include 'COMMON.TIME1'
3431 include 'COMMON.THREAD'
3432 include 'COMMON.SBRIDGE'
3433 include 'COMMON.CONTROL'
3434 include 'COMMON.MCM'
3435 include 'COMMON.MAP'
3436 include 'COMMON.HEADER'
3437 include 'COMMON.CSA'
3438 include 'COMMON.CHAIN'
3439 include 'COMMON.MUCA'
3441 include 'COMMON.FFIELD'
3442 include 'COMMON.SETUP'
3444 double precision seed,ran_number
3445 iseed=-dint(dabs(seed))
3446 if (iseed.eq.0) then
3447 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3448 & 'Random seed undefined. The program will stop.'
3449 write (*,'(/80(1h*)/20x,a/80(1h*))')
3450 & 'Random seed undefined. The program will stop.'
3452 call mpi_finalize(mpi_comm_world,ierr)
3454 stop 'Bad random seed.'
3457 if (fg_rank.eq.0) then
3461 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3462 OKRandom = prng_restart(me,iseed)
3465 tmp=65536.0d0**(4-i)
3466 iseed_array(i) = dint(seed/tmp)
3467 seed=seed-iseed_array(i)*tmp
3470 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3471 & (iseed_array(i),i=1,4)
3472 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3473 & (iseed_array(i),i=1,4)
3474 OKRandom = prng_restart(me,iseed_array)
3477 c r1 = prng_next(me)
3478 r1=ran_number(0.0D0,1.0D0)
3480 & write (iout,*) 'ran_num',r1
3481 if (r1.lt.0.0d0) OKRandom=.false.
3483 if (.not.OKRandom) then
3484 write (iout,*) 'PRNG IS NOT WORKING!!!'
3485 print *,'PRNG IS NOT WORKING!!!'
3488 call mpi_abort(mpi_comm_world,error_msg,ierr)
3491 write (iout,*) 'too many processors for parallel prng'
3492 write (*,*) 'too many processors for parallel prng'
3500 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3504 c----------------------------------------------------------------------
3505 subroutine read_klapaucjusz
3507 include 'DIMENSIONS'
3511 include 'COMMON.SETUP'
3512 include 'COMMON.CONTROL'
3513 include 'COMMON.HOMOLOGY'
3514 include 'COMMON.CHAIN'
3515 include 'COMMON.IOUNITS'
3517 include 'COMMON.GEO'
3518 include 'COMMON.INTERACT'
3519 include 'COMMON.NAMES'
3520 character*256 fragfile
3521 integer ninclust(maxclust),inclust(max_template,maxclust),
3522 & nresclust(maxclust),iresclust(maxres,maxclust),nclust
3525 character*24 model_ki_dist, model_ki_angle
3526 character*500 controlcard
3527 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
3528 & ik,ll,ii,kk,iistart,iishift,lim_xx
3529 double precision distal
3530 logical lprn /.true./
3536 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3537 double precision, dimension (max_template,maxres) :: rescore
3538 double precision, dimension (max_template,maxres) :: rescore2
3539 character*24 pdbfile,tpl_k_rescore
3542 c For new homol impl
3544 include 'COMMON.VAR'
3546 call getenv("FRAGFILE",fragfile)
3547 open(ientin,file=fragfile,status="old",err=10)
3548 read(ientin,*) constr_homology,nclust
3554 do k=1,constr_homology
3555 read(ientin,'(a)') pdbfile
3556 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3557 & pdbfile(:ilen(pdbfile))
3558 open(ipdbin,file=pdbfile,status='old',err=33)
3560 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3561 & pdbfile(:ilen(pdbfile))
3565 call readpdb_template(k)
3573 read(ientin,*) ninclust(i),nresclust(i)
3574 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3575 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3578 c Loop over clusters
3581 do ll = 1,ninclust(l)
3589 idomain(k,iresclust(i,l)+1) = 1
3591 idomain(k,iresclust(i,l)) = 1
3595 c Distance restraints
3598 C Copy the coordinates from reference coordinates (?)
3602 c write (iout,*) "c(",j,i,") =",c(j,i)
3605 call int_from_cart(.true.,.false.)
3606 call sc_loc_geom(.true.)
3608 thetaref(i)=theta(i)
3611 if (waga_dist.ne.0.0d0) then
3619 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3620 c write (iout,*) k,i,j,distal,dist2_cut
3622 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3623 & .and. distal.le.dist2_cut ) then
3629 c write (iout,*) "k",k
3630 c write (iout,*) "i",i," j",j," constr_homology",
3635 if (read2sigma) then
3638 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3640 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3641 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3642 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3644 if (odl(k,ii).le.dist_cut) then
3645 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3648 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3649 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3651 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3652 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3656 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3659 c l_homo(k,ii)=.false.
3666 c Theta, dihedral and SC retraints
3668 if (waga_angle.gt.0.0d0) then
3670 if (idomain(k,i).eq.0) then
3671 c sigma_dih(k,i)=0.0
3675 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3676 & rescore(k,i-2)+rescore(k,i-3))/4.0
3677 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3678 c & " sigma_dihed",sigma_dih(k,i)
3679 if (sigma_dih(k,i).ne.0)
3680 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3685 if (waga_theta.gt.0.0d0) then
3687 if (idomain(k,i).eq.0) then
3688 c sigma_theta(k,i)=0.0
3691 thetatpl(k,i)=thetaref(i)
3692 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3693 & rescore(k,i-2))/3.0
3694 if (sigma_theta(k,i).ne.0)
3695 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3699 if (waga_d.gt.0.0d0) then
3701 if (itype(i).eq.10) cycle
3702 if (idomain(k,i).eq.0 ) then
3709 sigma_d(k,i)=rescore(k,i)
3710 if (sigma_d(k,i).ne.0)
3711 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3712 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3718 c remove distance restraints not used in any model from the list
3719 c shift data in all arrays
3721 if (waga_dist.ne.0.0d0) then
3727 if (ii_in_use(ii).eq.0.and.liiflag) then
3731 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3732 & .not.liiflag.and.ii.eq.lim_odl) then
3733 if (ii.eq.lim_odl) then
3734 iishift=ii-iistart+1
3739 do ki=iistart,lim_odl-iishift
3740 ires_homo(ki)=ires_homo(ki+iishift)
3741 jres_homo(ki)=jres_homo(ki+iishift)
3742 ii_in_use(ki)=ii_in_use(ki+iishift)
3743 do k=1,constr_homology
3744 odl(k,ki)=odl(k,ki+iishift)
3745 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3746 l_homo(k,ki)=l_homo(k,ki+iishift)
3750 lim_odl=lim_odl-iishift
3757 10 stop "Error in fragment file"