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 readi(controlcard,"INTER_LIST_UPDATE",imatupdate,100)
175 call reada(controlcard,"T_BATH",t_bath,300.0d0)
176 minim=(index(controlcard,'MINIMIZE').gt.0)
177 dccart=(index(controlcard,'CART').gt.0)
178 overlapsc=(index(controlcard,'OVERLAP').gt.0)
179 overlapsc=.not.overlapsc
180 searchsc=(index(controlcard,'SEARCHSC').gt.0 .and.
181 & index(controlcard,'NOSEARCHSC').eq.0)
182 sideadd=(index(controlcard,'SIDEADD').gt.0)
183 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
184 mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
185 outpdb=(index(controlcard,'PDBOUT').gt.0)
186 outmol2=(index(controlcard,'MOL2OUT').gt.0)
187 pdbref=(index(controlcard,'PDBREF').gt.0)
188 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
189 indpdb=index(controlcard,'PDBSTART')
190 AFMlog=(index(controlcard,'AFM'))
191 selfguide=(index(controlcard,'SELFGUIDE'))
192 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
193 call readi(controlcard,'TUBEMOD',tubelog,0)
194 c write (iout,*) TUBElog,"TUBEMODE"
195 call readi(controlcard,'IPRINT',iprint,0)
196 C SHIELD keyword sets if the shielding effect of side-chains is used
197 C 0 denots no shielding is used all peptide are equally despite the
198 C solvent accesible area
199 C 1 the newly introduced function
200 C 2 reseved for further possible developement
201 call readi(controlcard,'SHIELD',shield_mode,0)
202 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
203 write(iout,*) "shield_mode",shield_mode
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 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
299 extconf=(index(controlcard,'EXTCONF').gt.0)
300 if (start_from_model) then
304 i2ndstr=index(controlcard,'USE_SEC_PRED')
305 gradout=index(controlcard,'GRADOUT').gt.0
306 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
307 C DISTCHAINMAX become obsolete for periodic boundry condition
308 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
309 C Reading the dimensions of box in x,y,z coordinates
310 call reada(controlcard,'BOXX',boxxsize,100.0d0)
311 call reada(controlcard,'BOXY',boxysize,100.0d0)
312 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
313 write(iout,*) "Periodic box dimensions",boxxsize,boxysize,boxzsize
314 c Cutoff range for interactions
315 call reada(controlcard,"R_CUT_INT",r_cut_int,25.0d0)
316 call reada(controlcard,"R_CUT_RESPA",r_cut_respa,2.0d0)
317 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
318 write (iout,*) "Cutoff on interactions",r_cut_int
320 & "Cutoff in switching short and long range interactions in RESPA",
322 write (iout,*) "lambda in switch function",rlamb
323 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
324 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
325 if (lipthick.gt.0.0d0) then
326 bordliptop=(boxzsize+lipthick)/2.0
327 bordlipbot=bordliptop-lipthick
329 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
330 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
331 buflipbot=bordlipbot+lipbufthick
332 bufliptop=bordliptop-lipbufthick
333 if ((lipbufthick*2.0d0).gt.lipthick)
334 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
336 write(iout,*) "bordliptop=",bordliptop
337 write(iout,*) "bordlipbot=",bordlipbot
338 write(iout,*) "bufliptop=",bufliptop
339 write(iout,*) "buflipbot=",buflipbot
340 write (iout,*) "SHIELD MODE",shield_mode
341 if (TUBElog.gt.0) then
342 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
343 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
344 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
345 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
346 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
347 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
348 buftubebot=bordtubebot+tubebufthick
349 buftubetop=bordtubetop-tubebufthick
351 if (me.eq.king .or. .not.out1file )
352 & write (iout,*) "DISTCHAINMAX",distchainmax
354 if(me.eq.king.or..not.out1file)
355 & write (iout,'(2a)') diagmeth(kdiag),
356 & ' routine used to diagonalize matrices.'
359 c--------------------------------------------------------------------------
360 subroutine read_REMDpar
366 include 'COMMON.IOUNITS'
367 include 'COMMON.TIME1'
370 include 'COMMON.LANGEVIN'
373 include 'COMMON.LANGEVIN.lang0.5diag'
375 include 'COMMON.LANGEVIN.lang0'
378 include 'COMMON.INTERACT'
379 include 'COMMON.NAMES'
381 include 'COMMON.REMD'
382 include 'COMMON.CONTROL'
383 include 'COMMON.SETUP'
385 character*320 controlcard
386 character*3200 controlcard1
387 integer iremd_m_total,i
389 if(me.eq.king.or..not.out1file)
390 & write (iout,*) "REMD setup"
392 call card_concat(controlcard)
393 call readi(controlcard,"NREP",nrep,3)
394 call readi(controlcard,"NSTEX",nstex,1000)
395 call reada(controlcard,"RETMIN",retmin,10.0d0)
396 call reada(controlcard,"RETMAX",retmax,1000.0d0)
397 mremdsync=(index(controlcard,'SYNC').gt.0)
398 call readi(controlcard,"NSYN",i_sync_step,100)
399 restart1file=(index(controlcard,'REST1FILE').gt.0)
400 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
401 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
402 if(max_cache_traj_use.gt.max_cache_traj)
403 & max_cache_traj_use=max_cache_traj
404 if(me.eq.king.or..not.out1file) then
405 cd if (traj1file) then
406 crc caching is in testing - NTWX is not ignored
407 cd write (iout,*) "NTWX value is ignored"
408 cd write (iout,*) " trajectory is stored to one file by master"
409 cd write (iout,*) " before exchange at NSTEX intervals"
411 write (iout,*) "NREP= ",nrep
412 write (iout,*) "NSTEX= ",nstex
413 write (iout,*) "SYNC= ",mremdsync
414 write (iout,*) "NSYN= ",i_sync_step
415 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
418 if (index(controlcard,'TLIST').gt.0) then
420 call card_concat(controlcard1)
421 read(controlcard1,*) (remd_t(i),i=1,nrep)
422 if(me.eq.king.or..not.out1file)
423 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
426 if (index(controlcard,'MLIST').gt.0) then
428 call card_concat(controlcard1)
429 read(controlcard1,*) (remd_m(i),i=1,nrep)
430 if(me.eq.king.or..not.out1file) then
431 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
434 iremd_m_total=iremd_m_total+remd_m(i)
436 write (iout,*) 'Total number of replicas ',iremd_m_total
441 & "Adaptive (PMF-biased) umbrella sampling will be run"
444 if(me.eq.king.or..not.out1file)
445 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
448 c--------------------------------------------------------------------------
449 subroutine read_MDpar
455 include 'COMMON.IOUNITS'
456 include 'COMMON.TIME1'
458 include 'COMMON.QRESTR'
460 include 'COMMON.LANGEVIN'
463 include 'COMMON.LANGEVIN.lang0.5diag'
465 include 'COMMON.LANGEVIN.lang0'
468 include 'COMMON.INTERACT'
469 include 'COMMON.NAMES'
471 include 'COMMON.SETUP'
472 include 'COMMON.CONTROL'
473 include 'COMMON.SPLITELE'
474 include 'COMMON.FFIELD'
476 character*320 controlcard
480 call card_concat(controlcard)
481 call readi(controlcard,"NSTEP",n_timestep,1000000)
482 call readi(controlcard,"NTWE",ntwe,100)
483 call readi(controlcard,"NTWX",ntwx,1000)
484 call reada(controlcard,"DT",d_time,1.0d-1)
485 call reada(controlcard,"DVMAX",dvmax,2.0d1)
486 call reada(controlcard,"DAMAX",damax,1.0d1)
487 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
488 call readi(controlcard,"LANG",lang,0)
489 RESPA = index(controlcard,"RESPA") .gt. 0
490 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
491 ntime_split0=ntime_split
492 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
493 ntime_split0=ntime_split
494 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
495 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
496 rest = index(controlcard,"REST").gt.0
497 tbf = index(controlcard,"TBF").gt.0
498 usampl = index(controlcard,"USAMPL").gt.0
499 scale_umb = index(controlcard,"SCALE_UMB").gt.0
500 adaptive = index(controlcard,"ADAPTIVE").gt.0
501 mdpdb = index(controlcard,"MDPDB").gt.0
502 call reada(controlcard,"T_BATH",t_bath,300.0d0)
503 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
504 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
505 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
506 if (count_reset_moment.eq.0) count_reset_moment=1000000000
507 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
508 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
509 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
510 if (count_reset_vel.eq.0) count_reset_vel=1000000000
511 large = index(controlcard,"LARGE").gt.0
512 print_compon = index(controlcard,"PRINT_COMPON").gt.0
513 rattle = index(controlcard,"RATTLE").gt.0
514 preminim = index(controlcard,"PREMINIM").gt.0
515 if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then
519 write (iout,*) "PREMINIM ",preminim
521 if (index(controlcard,'CART').gt.0) dccart=.true.
522 write (iout,*) "dccart ",dccart
523 write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR')
524 if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim
526 c if performing umbrella sampling, fragments constrained are read from the fragment file
529 write (iout,*) "Umbrella sampling will be run"
530 if (scale_umb.and.adaptive) then
531 write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
532 write (iout,*) "Select one of those and re-run the job."
535 if (scale_umb) write (iout,*)
536 &"Umbrella-restraint force constants will be scaled by temperature"
540 if(me.eq.king.or..not.out1file) then
542 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
544 write (iout,'(a)') "The units are:"
545 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
546 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
547 & " acceleration: angstrom/(48.9 fs)**2"
548 write (iout,'(a)') "energy: kcal/mol, temperature: K"
550 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
551 write (iout,'(a60,f10.5,a)')
552 & "Initial time step of numerical integration:",d_time,
554 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
555 write (iout,'(a60,f10.5,a)') "Cutoff on interactions",r_cut_int,
557 write(iout,'(a60,i5)')"Frequency of updating interaction list",
560 write (iout,'(2a,i4,a)')
561 & "A-MTS algorithm used; initial time step for fast-varying",
562 & " short-range forces split into",ntime_split," steps."
563 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
564 & r_cut_respa," lambda",rlamb
566 write (iout,'(2a,f10.5)')
567 & "Maximum acceleration threshold to reduce the time step",
568 & "/increase split number:",damax
569 write (iout,'(2a,f10.5)')
570 & "Maximum predicted energy drift to reduce the timestep",
571 & "/increase split number:",edriftmax
572 write (iout,'(a60,f10.5)')
573 & "Maximum velocity threshold to reduce velocities:",dvmax
574 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
575 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
576 if (rattle) write (iout,'(a60)')
577 & "Rattle algorithm used to constrain the virtual bonds"
581 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
582 call reada(controlcard,"RWAT",rwat,1.4d0)
583 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
584 surfarea=index(controlcard,"SURFAREA").gt.0
585 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
586 if(me.eq.king.or..not.out1file)then
587 write (iout,'(/a,$)') "Langevin dynamics calculation"
590 & " with direct integration of Langevin equations"
591 else if (lang.eq.2) then
592 write (iout,'(a/)') " with TINKER stochasic MD integrator"
593 else if (lang.eq.3) then
594 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
595 else if (lang.eq.4) then
596 write (iout,'(a/)') " in overdamped mode"
598 write (iout,'(//a,i5)')
599 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
602 write (iout,'(a60,f10.5)') "Temperature:",t_bath
603 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
604 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
605 write (iout,'(a60,f10.5)')
606 & "Scaling factor of the friction forces:",scal_fric
607 if (surfarea) write (iout,'(2a,i10,a)')
608 & "Friction coefficients will be scaled by solvent-accessible",
609 & " surface area every",reset_fricmat," steps."
611 c Calculate friction coefficients and bounds of stochastic forces
612 eta=6*pi*cPoise*etawat
613 if(me.eq.king.or..not.out1file)
614 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
616 gamp=scal_fric*(pstok+rwat)*eta
617 stdfp=dsqrt(2*Rb*t_bath/d_time)
619 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
620 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
622 if(me.eq.king.or..not.out1file)then
623 write (iout,'(/2a/)')
624 & "Radii of site types and friction coefficients and std's of",
625 & " stochastic forces of fully exposed sites"
626 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
628 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
629 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
633 if(me.eq.king.or..not.out1file)then
634 write (iout,'(a)') "Berendsen bath calculation"
635 write (iout,'(a60,f10.5)') "Temperature:",t_bath
636 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
638 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
639 & count_reset_moment," steps"
641 & write (iout,'(a,i10,a)')
642 & "Velocities will be reset at random every",count_reset_vel,
646 if(me.eq.king.or..not.out1file)
647 & write (iout,'(a31)') "Microcanonical mode calculation"
649 if(me.eq.king.or..not.out1file)then
650 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
652 write(iout,*) "MD running with constraints."
653 write(iout,*) "Equilibration time ", eq_time, " mtus."
654 write(iout,*) "Constraining ", nfrag," fragments."
655 write(iout,*) "Length of each fragment, weight and q0:"
657 write (iout,*) "Set of restraints #",iset
659 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
660 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
662 write(iout,*) "constraints between ", npair, "fragments."
663 write(iout,*) "constraint pairs, weights and q0:"
665 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
666 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
668 write(iout,*) "angle constraints within ", nfrag_back,
669 & "backbone fragments."
671 write(iout,*) "fragment, weights, q0:"
673 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
674 & ifrag_back(2,i,iset),
675 & wfrag_back(1,i,iset),qin_back(1,i,iset),
676 & wfrag_back(2,i,iset),qin_back(2,i,iset),
677 & wfrag_back(3,i,iset),qin_back(3,i,iset)
680 write(iout,*) "fragment, weights:"
682 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
683 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
684 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
688 iset=mod(kolor,nset)+1
691 if(me.eq.king.or..not.out1file)
692 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
695 c------------------------------------------------------------------------------
698 C Read molecular data.
704 integer error_msg,ierror,ierr,ierrcode
706 include 'COMMON.IOUNITS'
709 include 'COMMON.INTERACT'
710 include 'COMMON.LOCAL'
711 include 'COMMON.NAMES'
712 include 'COMMON.CHAIN'
713 include 'COMMON.FFIELD'
714 include 'COMMON.SBRIDGE'
715 include 'COMMON.HEADER'
716 include 'COMMON.CONTROL'
717 include 'COMMON.SAXS'
718 include 'COMMON.DBASE'
719 include 'COMMON.THREAD'
720 include 'COMMON.CONTACTS'
721 include 'COMMON.TORCNSTR'
722 include 'COMMON.TIME1'
723 include 'COMMON.BOUNDS'
725 include 'COMMON.SETUP'
726 include 'COMMON.SHIELD'
727 character*4 sequence(maxres)
729 double precision x(maxvar)
730 character*256 pdbfile
731 character*400 weightcard
732 character*80 weightcard_t,ucase
733 integer itype_pdb(maxres)
734 common /pizda/ itype_pdb
735 logical seq_comp,fail
736 double precision energia(0:n_ene)
737 double precision secprob(3,maxdih_constr)
739 double precision phihel,phibet,sigmahel,sigmabet
740 integer iti,nsi,maxsi
744 integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp,itemp
745 double precision sumv
747 C Read PDB structure if applicable
749 if (indpdb.gt.0 .or. pdbref) then
750 read(inp,'(a)') pdbfile
751 if(me.eq.king.or..not.out1file)
752 & write (iout,'(2a)') 'PDB data will be read from file ',
753 & pdbfile(:ilen(pdbfile))
754 open(ipdbin,file=pdbfile,status='old',err=33)
756 33 write (iout,'(a)') 'Error opening PDB file.'
767 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
768 & (crefjlee(j,i+nres),j=1,3)
771 if(me.eq.king.or..not.out1file)
772 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
773 & ' nstart_sup=',nstart_sup
775 itype_pdb(i)=itype(i)
779 nct=nstart_sup+nsup-1
780 call contact(.false.,ncont_ref,icont_ref,co)
783 if(me.eq.king.or..not.out1file)
784 & write(iout,*)'Adding sidechains'
788 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
791 do while (fail.and.nsi.le.maxsi)
792 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
795 if(fail) write(iout,*)'Adding sidechain failed for res ',
796 & i,' after ',nsi,' trials'
801 if (indpdb.eq.0) then
802 C Read sequence if not taken from the pdb file.
804 c print *,'nres=',nres
805 if (iscode.gt.0) then
806 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
808 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
810 C Convert sequence to numeric code
812 itype(i)=rescode(i,sequence(i),iscode)
816 c print '(20i4)',(itype(i),i=1,nres)
819 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
821 if (itype(i).eq.ntyp1) then
825 else if (iabs(itype(i+1)).ne.20) then
827 else if (iabs(itype(i)).ne.20) then
834 if(me.eq.king.or..not.out1file)then
835 write (iout,*) "ITEL"
837 write (iout,*) i,itype(i),itel(i)
839 c print *,'Call Read_Bridge.'
843 cd print *,'NNT=',NNT,' NCT=',NCT
844 call seq2chains(nres,itype,nchain,chain_length,chain_border,
847 chain_border1(2,1)=chain_border(2,1)+1
849 chain_border1(1,i)=chain_border(1,i)-1
850 chain_border1(2,i)=chain_border(2,i)+1
852 if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
853 chain_border1(2,nchain)=nres
854 write(iout,*) "nres",nres," nchain",nchain
856 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
857 & chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
859 call chain_symmetry(nchain,nres,itype,chain_border,
860 & chain_length,npermchain,tabpermchain)
862 c write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
865 write(iout,*) "residue permutations"
867 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
870 if (itype(1).eq.ntyp1) nnt=2
871 if (itype(nres).eq.ntyp1) nct=nct-1
872 write (iout,*) "nnt",nnt," nct",nct
875 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
876 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
878 c print*, 'init_dfa_vars finished!'
880 c print*, 'read_dfa_info finished!'
884 if(me.eq.king.or..not.out1file)
885 & write (iout,'(a,i3)') 'nsup=',nsup
887 if (nsup.le.(nct-nnt+1)) then
888 do i=0,nct-nnt+1-nsup
889 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
895 & 'Error - sequences to be superposed do not match.'
898 do i=0,nsup-(nct-nnt+1)
899 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
901 nstart_sup=nstart_sup+i
907 & 'Error - sequences to be superposed do not match.'
910 if (nsup.eq.0) nsup=nct-nnt
911 if (nstart_sup.eq.0) nstart_sup=nnt
912 if (nstart_seq.eq.0) nstart_seq=nnt
913 if(me.eq.king.or..not.out1file)
914 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
915 & ' nstart_seq=',nstart_seq
918 C 8/13/98 Set limits to generating the dihedral angles
923 read (inp,*) ndih_constr
924 write (iout,*) "ndih_constr",ndih_constr
925 if (ndih_constr.gt.0) then
928 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
931 if(me.eq.king.or..not.out1file)then
934 & 'There are',ndih_constr,' restraints on gamma angles.'
936 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
943 phi0(i)=deg2rad*phi0(i)
944 drange(i)=deg2rad*drange(i)
946 C if(me.eq.king.or..not.out1file)
947 C & write (iout,*) 'FTORS',ftors
950 phibound(1,ii) = phi0(i)-drange(i)
951 phibound(2,ii) = phi0(i)+drange(i)
954 if (me.eq.king .or. .not.out1file) then
956 write (iout,'(a)') 'Boundaries in gamma angle sampling:'
958 write (iout,'(a3,i5,2f10.1)')
959 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
964 else if (ndih_constr.lt.0) then
966 call card_concat(weightcard)
967 call reada(weightcard,"PHIHEL",phihel,50.0D0)
968 call reada(weightcard,"PHIBET",phibet,180.0D0)
969 call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
970 call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
971 call reada(weightcard,"WDIHC",wdihc,0.591D0)
972 write (iout,*) "Weight of dihedral angle restraints",wdihc
973 read(inp,'(9x,3f7.3)')
974 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
975 write (iout,*) "The secprob array"
977 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
981 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
982 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
983 ndih_constr=ndih_constr+1
984 idih_constr(ndih_constr)=i
987 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
988 sumv=sumv+vpsipred(j,ndih_constr)
991 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
993 phibound(1,ndih_constr)=phihel*deg2rad
994 phibound(2,ndih_constr)=phibet*deg2rad
995 sdihed(1,ndih_constr)=sigmahel*deg2rad
996 sdihed(2,ndih_constr)=sigmabet*deg2rad
1000 if(me.eq.king.or..not.out1file)then
1003 & 'There are',ndih_constr,
1004 & ' bimodal restraints on gamma angles.'
1006 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
1007 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
1008 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
1009 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
1010 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
1011 & (vpsipred(j,i),j=1,3)
1018 C first setting the theta boundaries to 0 to pi
1019 C this mean that there is no energy penalty for any angle occuring this can be applied
1020 C for generate random conformation but is not implemented in this way
1023 C thetabound(2,i)=pi
1025 C begin reading theta constrains this is quartic constrains allowing to
1026 C have smooth second derivative
1027 if (with_theta_constr) then
1028 C with_theta_constr is keyword allowing for occurance of theta constrains
1029 read (inp,*) ntheta_constr
1030 C ntheta_constr is the number of theta constrains
1031 if (ntheta_constr.gt.0) then
1032 C read (inp,*) ftors
1033 read (inp,*) (itheta_constr(i),theta_constr0(i),
1034 & theta_drange(i),for_thet_constr(i),
1035 & i=1,ntheta_constr)
1036 C the above code reads from 1 to ntheta_constr
1037 C itheta_constr(i) residue i for which is theta_constr
1038 C theta_constr0 the global minimum value
1039 C theta_drange is range for which there is no energy penalty
1040 C for_thet_constr is the force constant for quartic energy penalty
1042 if(me.eq.king.or..not.out1file)then
1044 & 'There are',ntheta_constr,' constraints on phi angles.'
1045 do i=1,ntheta_constr
1046 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1048 & for_thet_constr(i)
1051 do i=1,ntheta_constr
1052 theta_constr0(i)=deg2rad*theta_constr0(i)
1053 theta_drange(i)=deg2rad*theta_drange(i)
1055 C if(me.eq.king.or..not.out1file)
1056 C & write (iout,*) 'FTORS',ftors
1057 C do i=1,ntheta_constr
1058 C ii = itheta_constr(i)
1059 C thetabound(1,ii) = phi0(i)-drange(i)
1060 C thetabound(2,ii) = phi0(i)+drange(i)
1062 endif ! ntheta_constr.gt.0
1063 endif! with_theta_constr
1065 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1066 C write (iout,*) "with_dihed_constr ",with_dihed_constr
1067 c--- Zscore rms -------
1068 if (nz_start.eq.0) nz_start=nnt
1069 if (nz_end.eq.0 .and. nsup.gt.0) then
1071 else if (nz_end.eq.0) then
1074 if(me.eq.king.or..not.out1file)then
1075 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1076 write (iout,*) 'IZ_SC=',iz_sc
1078 c----------------------
1081 if (.not.pdbref) then
1082 call read_angles(inp,*38)
1085 38 write (iout,'(a)') 'Error reading reference structure.'
1087 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1088 stop 'Error reading reference structure'
1090 39 call chainbuild_extconf
1092 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1101 call contact(.true.,ncont_ref,icont_ref,co)
1105 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1107 if (constr_dist.gt.0) call read_dist_constr
1108 c write (iout,*) "After read_dist_constr nhpb",nhpb
1109 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1111 call NMRpeak_partition
1112 if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1114 if(me.eq.king.or..not.out1file)
1115 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1118 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1120 if(me.eq.king.or..not.out1file)
1121 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1122 & icont_ref(1,i),' ',
1123 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1126 write (iout,*) "calling read_saxs_consrtr",nsaxs
1127 if (nsaxs.gt.0) call read_saxs_constr
1129 if (constr_homology.gt.0) then
1130 call read_constr_homology
1131 if (indpdb.gt.0 .or. pdbref) then
1134 c(j,i)=crefjlee(j,i)
1135 cref(j,i)=crefjlee(j,i)
1140 write (iout,*) "sc_loc_geom: Array C"
1142 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1143 & (c(j,i+nres),j=1,3)
1145 write (iout,*) "Array Cref"
1147 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1148 & (cref(j,i+nres),j=1,3)
1151 call int_from_cart1(.false.)
1152 call sc_loc_geom(.false.)
1154 thetaref(i)=theta(i)
1159 dc(j,i)=c(j,i+1)-c(j,i)
1160 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1165 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1166 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1171 if (start_from_model) then
1174 read(inp,'(a)',end=332,err=332) pdbfile
1175 if (me.eq.king .or. .not. out1file)
1176 & write (iout,'(a,5x,a)') 'Opening PDB file',
1177 & pdbfile(:ilen(pdbfile))
1178 open(ipdbin,file=pdbfile,status='old',err=336)
1180 336 write (iout,'(a,5x,a)') 'Error opening PDB file',
1181 & pdbfile(:ilen(pdbfile))
1189 if (nres.ge.nres_temp) then
1190 nmodel_start=nmodel_start+1
1191 pdbfiles_chomo(nmodel_start)=pdbfile
1194 chomo(j,i,nmodel_start)=c(j,i)
1200 c call gen_rand_conf(itemp,*115)
1201 c nmodel_start=nmodel_start+1
1204 c chomo(j,i,nmodel_start)=c(j,i)
1208 115 if (me.eq.king .or. .not. out1file)
1209 & write (iout,'(a,2i5,1x,a)')
1210 & "Different number of residues",nres_temp,nres,
1217 if (nmodel_start.eq.0) then
1218 if (me.eq.king .or. .not. out1file)
1219 & write (iout,'(a)')
1220 & "No valid starting model found START_FROM_MODELS is OFF"
1221 start_from_model=.false.
1223 write (iout,*) "nmodel_start",nmodel_start
1229 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1230 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1231 & modecalc.ne.10) then
1232 C If input structure hasn't been supplied from the PDB file read or generate
1234 if (iranconf.eq.0 .and. .not. extconf .and. .not.
1235 & start_from_model) then
1236 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1237 & write (iout,'(a)') 'Initial geometry will be read in.'
1239 read(inp,'(8f10.5)',end=36,err=36)
1240 & ((c(l,k),l=1,3),k=1,nres),
1241 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1242 if (nnt.gt.1) c(:,nres+1)=c(:,1)
1243 if (nct.lt.nres) c(:,2*nres)=c(:,nres)
1244 c write (iout,*) "Exit READ_CART"
1245 c write (iout,'(8f10.5)')
1246 c & ((c(l,k),l=1,3),k=1,nres),
1247 c & ((c(l,k+nres),l=1,3),k=nnt,nct)
1254 dc(j,i)=c(j,i+1)-c(j,i)
1255 c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1259 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1261 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1262 c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1267 c dc_norm(j,i+nres)=0.0d0
1271 call int_from_cart1(.true.)
1272 write (iout,*) "Finish INT_TO_CART"
1273 c write (iout,*) "DC"
1275 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1276 c & (dc(j,i+nres),j=1,3)
1282 call read_angles(inp,*36)
1284 call chainbuild_extconf
1287 36 write (iout,'(a)') 'Error reading angle file.'
1289 call mpi_finalize( MPI_COMM_WORLD,IERR )
1291 stop 'Error reading angle file.'
1293 else if (extconf) then
1294 if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1295 & write (iout,'(a)') 'Extended chain initial geometry.'
1297 theta(i)=90d0*deg2rad
1300 phi(i)=180d0*deg2rad
1303 alph(i)=110d0*deg2rad
1306 omeg(i)=-120d0*deg2rad
1307 if (itype(i).le.0) omeg(i)=-omeg(i)
1310 call chainbuild_extconf
1311 else if (.not. start_from_model) then
1312 if(me.eq.king.or..not.out1file)
1313 & write (iout,'(a)') 'Random-generated initial geometry.'
1316 if (me.eq.king .or. fg_rank.eq.0 .and. (
1317 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1321 call gen_rand_conf(itmp,*30)
1323 30 write (iout,*) 'Failed to generate random conformation',
1324 & ', itrial=',itrial
1325 write (*,*) 'Processor:',me,
1326 & ' Failed to generate random conformation',
1336 write (iout,'(a,i3,a)') 'Processor:',me,
1337 & ' error in generating random conformation.'
1338 write (*,'(a,i3,a)') 'Processor:',me,
1339 & ' error in generating random conformation.'
1342 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1347 & ' error in generating random conformation.'
1352 elseif (modecalc.eq.4) then
1353 read (inp,'(a)') intinname
1354 open (intin,file=intinname,status='old',err=333)
1355 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1356 & write (iout,'(a)') 'intinname',intinname
1357 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1359 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1361 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1363 stop 'Error opening angle file.'
1367 C Generate distance constraints, if the PDB structure is to be regularized.
1368 if (nthread.gt.0) then
1369 call read_threadbase
1372 if (me.eq.king .or. .not. out1file)
1374 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1375 write (iout,'(/a,i3,a)')
1376 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1377 write (iout,'(20i4)') (iss(i),i=1,ns)
1379 write(iout,*)"Running with dynamic disulfide-bond formation"
1381 write (iout,'(/a/)') 'Pre-formed links are:'
1387 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1388 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1394 if (ns.gt.0.and.dyn_ss) then
1398 forcon(i-nss)=forcon(i)
1405 dyn_ss_mask(iss(i))=.true.
1410 c write (iout,*) "DC"
1412 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1413 c & (dc(j,i+nres),j=1,3)
1415 if (i2ndstr.gt.0) call secstrp2dihc
1416 c call geom_to_var(nvar,x)
1417 c call etotal(energia(0))
1418 c call enerprint(energia(0))
1419 c call briefout(0,etot)
1421 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1422 cd write (iout,'(a)') 'Variable list:'
1423 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1425 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1426 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1427 & 'Processor',myrank,': end reading molecular data.'
1432 c--------------------------------------------------------------------------
1433 logical function seq_comp(itypea,itypeb,length)
1435 integer length,itypea(length),itypeb(length)
1438 if (itypea(i).ne.itypeb(i)) then
1446 c-----------------------------------------------------------------------------
1447 subroutine read_bridge
1448 C Read information about disulfide bridges.
1450 include 'DIMENSIONS'
1455 include 'COMMON.IOUNITS'
1456 include 'COMMON.GEO'
1457 include 'COMMON.VAR'
1458 include 'COMMON.INTERACT'
1459 include 'COMMON.LOCAL'
1460 include 'COMMON.NAMES'
1461 include 'COMMON.CHAIN'
1462 include 'COMMON.FFIELD'
1463 include 'COMMON.SBRIDGE'
1464 include 'COMMON.HEADER'
1465 include 'COMMON.CONTROL'
1466 include 'COMMON.DBASE'
1467 include 'COMMON.THREAD'
1468 include 'COMMON.TIME1'
1469 include 'COMMON.SETUP'
1471 C Read bridging residues.
1472 read (inp,*) ns,(iss(i),i=1,ns)
1473 c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
1480 if(me.eq.king.or..not.out1file)
1481 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1482 C Check whether the specified bridging residues are cystines.
1484 if (itype(iss(i)).ne.1) then
1485 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1486 & 'Do you REALLY think that the residue ',
1487 & restyp(itype(iss(i))),i,
1488 & ' can form a disulfide bridge?!!!'
1489 write (*,'(2a,i3,a)')
1490 & 'Do you REALLY think that the residue ',
1491 & restyp(itype(iss(i))),i,
1492 & ' can form a disulfide bridge?!!!'
1494 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1499 C Read preformed bridges.
1501 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1503 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1506 C Check if the residues involved in bridges are in the specified list of
1507 C bridging residues.
1510 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1511 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1512 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1513 & ' contains residues present in other pairs.'
1514 write (*,'(a,i3,a)') 'Disulfide pair',i,
1515 & ' contains residues present in other pairs.'
1517 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1523 if (ihpb(i).eq.iss(j)) goto 10
1525 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1528 if (jhpb(i).eq.iss(j)) goto 20
1530 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1536 ihpb(i)=ihpb(i)+nres
1537 jhpb(i)=jhpb(i)+nres
1543 c----------------------------------------------------------------------------
1544 subroutine read_x(kanal,*)
1546 include 'DIMENSIONS'
1547 include 'COMMON.GEO'
1548 include 'COMMON.VAR'
1549 include 'COMMON.CHAIN'
1550 include 'COMMON.IOUNITS'
1551 include 'COMMON.CONTROL'
1552 include 'COMMON.LOCAL'
1553 include 'COMMON.INTERACT'
1554 integer i,j,k,l,kanal
1555 c Read coordinates from input
1557 read(kanal,'(8f10.5)',end=10,err=10)
1558 & ((c(l,k),l=1,3),k=1,nres),
1559 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1562 c(j,2*nres)=c(j,nres)
1564 call int_from_cart1(.false.)
1567 dc(j,i)=c(j,i+1)-c(j,i)
1568 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1572 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1574 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1575 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1583 c----------------------------------------------------------------------------
1584 subroutine read_threadbase
1586 include 'DIMENSIONS'
1587 include 'COMMON.IOUNITS'
1588 include 'COMMON.GEO'
1589 include 'COMMON.VAR'
1590 include 'COMMON.INTERACT'
1591 include 'COMMON.LOCAL'
1592 include 'COMMON.NAMES'
1593 include 'COMMON.CHAIN'
1594 include 'COMMON.FFIELD'
1595 include 'COMMON.SBRIDGE'
1596 include 'COMMON.HEADER'
1597 include 'COMMON.CONTROL'
1598 include 'COMMON.DBASE'
1599 include 'COMMON.THREAD'
1600 include 'COMMON.TIME1'
1602 double precision dist
1603 C Read pattern database for threading.
1604 read (icbase,*) nseq
1606 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1607 & nres_base(2,i),nres_base(3,i)
1608 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1610 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1611 c & nres_base(2,i),nres_base(3,i)
1612 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1616 if (weidis.eq.0.0D0) weidis=0.1D0
1625 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1626 write (iout,'(a,i5)') 'nexcl: ',nexcl
1627 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1630 c------------------------------------------------------------------------------
1631 subroutine setup_var
1633 include 'DIMENSIONS'
1634 include 'COMMON.IOUNITS'
1635 include 'COMMON.GEO'
1636 include 'COMMON.VAR'
1637 include 'COMMON.INTERACT'
1638 include 'COMMON.LOCAL'
1639 include 'COMMON.NAMES'
1640 include 'COMMON.CHAIN'
1641 include 'COMMON.FFIELD'
1642 include 'COMMON.SBRIDGE'
1643 include 'COMMON.HEADER'
1644 include 'COMMON.CONTROL'
1645 include 'COMMON.DBASE'
1646 include 'COMMON.THREAD'
1647 include 'COMMON.TIME1'
1649 C Set up variable list.
1655 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1657 ialph(i,1)=nvar+nside
1661 if (indphi.gt.0) then
1663 else if (indback.gt.0) then
1670 c----------------------------------------------------------------------------
1671 subroutine gen_dist_constr
1672 C Generate CA distance constraints.
1674 include 'DIMENSIONS'
1675 include 'COMMON.IOUNITS'
1676 include 'COMMON.GEO'
1677 include 'COMMON.VAR'
1678 include 'COMMON.INTERACT'
1679 include 'COMMON.LOCAL'
1680 include 'COMMON.NAMES'
1681 include 'COMMON.CHAIN'
1682 include 'COMMON.FFIELD'
1683 include 'COMMON.SBRIDGE'
1684 include 'COMMON.HEADER'
1685 include 'COMMON.CONTROL'
1686 include 'COMMON.DBASE'
1687 include 'COMMON.THREAD'
1688 include 'COMMON.SPLITELE'
1689 include 'COMMON.TIME1'
1690 integer i,j,itype_pdb(maxres)
1691 common /pizda/ itype_pdb
1693 double precision dist
1695 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1696 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1697 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1699 do i=nstart_sup,nstart_sup+nsup-1
1700 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1701 cd & ' seq_pdb', restyp(itype_pdb(i))
1702 do j=i+2,nstart_sup+nsup-1
1703 c 5/24/2020 Adam: Cutoff included to reduce array size
1705 if (dd.gt.r_cut_int) cycle
1707 ihpb(nhpb)=i+nstart_seq-nstart_sup
1708 jhpb(nhpb)=j+nstart_seq-nstart_sup
1713 cd write (iout,'(a)') 'Distance constraints:'
1718 cd if (ii.gt.nres) then
1723 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1724 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1725 cd & dhpb(i),forcon(i)
1729 c----------------------------------------------------------------------------
1732 include 'DIMENSIONS'
1733 include 'COMMON.MAP'
1734 include 'COMMON.IOUNITS'
1736 character*3 angid(4) /'THE','PHI','ALP','OME'/
1737 character*80 mapcard,ucase
1739 read (inp,'(a)') mapcard
1740 mapcard=ucase(mapcard)
1741 if (index(mapcard,'PHI').gt.0) then
1743 else if (index(mapcard,'THE').gt.0) then
1745 else if (index(mapcard,'ALP').gt.0) then
1747 else if (index(mapcard,'OME').gt.0) then
1750 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1751 stop 'Error - illegal variable spec in MAP card.'
1753 call readi (mapcard,'RES1',res1(imap),0)
1754 call readi (mapcard,'RES2',res2(imap),0)
1755 if (res1(imap).eq.0) then
1756 res1(imap)=res2(imap)
1757 else if (res2(imap).eq.0) then
1758 res2(imap)=res1(imap)
1760 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1762 & 'Error - illegal definition of variable group in MAP.'
1763 stop 'Error - illegal definition of variable group in MAP.'
1765 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1766 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1767 call readi(mapcard,'NSTEP',nstep(imap),0)
1768 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1770 & 'Illegal boundary and/or step size specification in MAP.'
1771 stop 'Illegal boundary and/or step size specification in MAP.'
1776 c----------------------------------------------------------------------------
1779 include 'DIMENSIONS'
1780 include 'COMMON.IOUNITS'
1781 include 'COMMON.GEO'
1782 include 'COMMON.CSA'
1783 include 'COMMON.BANK'
1784 include 'COMMON.CONTROL'
1786 character*620 mcmcard
1787 call card_concat(mcmcard)
1789 call readi(mcmcard,'NCONF',nconf,50)
1790 call readi(mcmcard,'NADD',nadd,0)
1791 call readi(mcmcard,'JSTART',jstart,1)
1792 call readi(mcmcard,'JEND',jend,1)
1793 call readi(mcmcard,'NSTMAX',nstmax,500000)
1794 call readi(mcmcard,'N0',n0,1)
1795 call readi(mcmcard,'N1',n1,6)
1796 call readi(mcmcard,'N2',n2,4)
1797 call readi(mcmcard,'N3',n3,0)
1798 call readi(mcmcard,'N4',n4,0)
1799 call readi(mcmcard,'N5',n5,0)
1800 call readi(mcmcard,'N6',n6,10)
1801 call readi(mcmcard,'N7',n7,0)
1802 call readi(mcmcard,'N8',n8,0)
1803 call readi(mcmcard,'N9',n9,0)
1804 call readi(mcmcard,'N14',n14,0)
1805 call readi(mcmcard,'N15',n15,0)
1806 call readi(mcmcard,'N16',n16,0)
1807 call readi(mcmcard,'N17',n17,0)
1808 call readi(mcmcard,'N18',n18,0)
1810 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1812 call readi(mcmcard,'NDIFF',ndiff,2)
1813 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1814 call readi(mcmcard,'IS1',is1,1)
1815 call readi(mcmcard,'IS2',is2,8)
1816 call readi(mcmcard,'NRAN0',nran0,4)
1817 call readi(mcmcard,'NRAN1',nran1,2)
1818 call readi(mcmcard,'IRR',irr,1)
1819 call readi(mcmcard,'NSEED',nseed,20)
1820 call readi(mcmcard,'NTOTAL',ntotal,10000)
1821 call reada(mcmcard,'CUT1',cut1,2.0d0)
1822 call reada(mcmcard,'CUT2',cut2,5.0d0)
1823 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1824 call readi(mcmcard,'ICMAX',icmax,3)
1825 call readi(mcmcard,'IRESTART',irestart,0)
1826 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1829 call reada(mcmcard,'DELE',dele,20.0d0)
1830 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1831 call readi(mcmcard,'IREF',iref,0)
1832 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1833 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1834 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1835 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1836 write (iout,*) "NCONF_IN",nconf_in
1839 c----------------------------------------------------------------------------
1842 include 'DIMENSIONS'
1843 include 'COMMON.MCM'
1844 include 'COMMON.MCE'
1845 include 'COMMON.IOUNITS'
1847 character*320 mcmcard
1849 call card_concat(mcmcard)
1850 call readi(mcmcard,'MAXACC',maxacc,100)
1851 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1852 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1853 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1854 call readi(mcmcard,'MAXREPM',maxrepm,200)
1855 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1856 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1857 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1858 call reada(mcmcard,'E_UP',e_up,5.0D0)
1859 call reada(mcmcard,'DELTE',delte,0.1D0)
1860 call readi(mcmcard,'NSWEEP',nsweep,5)
1861 call readi(mcmcard,'NSTEPH',nsteph,0)
1862 call readi(mcmcard,'NSTEPC',nstepc,0)
1863 call reada(mcmcard,'TMIN',tmin,298.0D0)
1864 call reada(mcmcard,'TMAX',tmax,298.0D0)
1865 call readi(mcmcard,'NWINDOW',nwindow,0)
1866 call readi(mcmcard,'PRINT_MC',print_mc,0)
1867 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1868 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1869 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1870 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1871 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1872 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1873 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1874 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1875 if (nwindow.gt.0) then
1876 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1878 winlen(i)=winend(i)-winstart(i)+1
1881 if (tmax.lt.tmin) tmax=tmin
1882 if (tmax.eq.tmin) then
1886 if (nstepc.gt.0 .and. nsteph.gt.0) then
1887 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1888 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1890 C Probabilities of different move types
1891 sumpro_type(0)=0.0D0
1892 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1893 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1894 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1895 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1896 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1897 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1898 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1900 print *,'i',i,' sumprotype',sumpro_type(i)
1901 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1902 print *,'i',i,' sumprotype',sumpro_type(i)
1906 c----------------------------------------------------------------------------
1907 subroutine read_minim
1909 include 'DIMENSIONS'
1910 include 'COMMON.MINIM'
1911 include 'COMMON.IOUNITS'
1913 character*320 minimcard
1914 call card_concat(minimcard)
1915 call readi(minimcard,'MAXMIN',maxmin,2000)
1916 call readi(minimcard,'MAXFUN',maxfun,5000)
1917 call readi(minimcard,'MINMIN',minmin,maxmin)
1918 call readi(minimcard,'MINFUN',minfun,maxmin)
1919 call reada(minimcard,'TOLF',tolf,1.0D-2)
1920 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1921 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1922 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1923 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1924 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1925 & 'Options in energy minimization:'
1926 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1927 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1928 & 'MinMin:',MinMin,' MinFun:',MinFun,
1929 & ' TolF:',TolF,' RTolF:',RTolF
1932 c----------------------------------------------------------------------------
1933 subroutine read_angles(kanal,*)
1935 include 'DIMENSIONS'
1936 include 'COMMON.GEO'
1937 include 'COMMON.VAR'
1938 include 'COMMON.CHAIN'
1939 include 'COMMON.IOUNITS'
1940 include 'COMMON.CONTROL'
1942 c Read angles from input
1944 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1945 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1946 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1947 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1950 c 9/7/01 avoid 180 deg valence angle
1951 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1953 theta(i)=deg2rad*theta(i)
1954 phi(i)=deg2rad*phi(i)
1955 alph(i)=deg2rad*alph(i)
1956 omeg(i)=deg2rad*omeg(i)
1961 c----------------------------------------------------------------------------
1962 subroutine reada(rekord,lancuch,wartosc,default)
1964 character*(*) rekord,lancuch
1965 double precision wartosc,default
1968 iread=index(rekord,lancuch)
1969 if (iread.eq.0) then
1973 iread=iread+ilen(lancuch)+1
1974 read (rekord(iread:),*,err=10,end=10) wartosc
1979 c----------------------------------------------------------------------------
1980 subroutine readi(rekord,lancuch,wartosc,default)
1982 character*(*) rekord,lancuch
1983 integer wartosc,default
1986 iread=index(rekord,lancuch)
1987 if (iread.eq.0) then
1991 iread=iread+ilen(lancuch)+1
1992 read (rekord(iread:),*,err=10,end=10) wartosc
1997 c----------------------------------------------------------------------------
1998 subroutine multreadi(rekord,lancuch,tablica,dim,default)
2001 integer tablica(dim),default
2002 character*(*) rekord,lancuch
2009 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2010 if (iread.eq.0) return
2011 iread=iread+ilen(lancuch)+1
2012 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2015 c----------------------------------------------------------------------------
2016 subroutine multreada(rekord,lancuch,tablica,dim,default)
2019 double precision tablica(dim),default
2020 character*(*) rekord,lancuch
2027 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2028 if (iread.eq.0) return
2029 iread=iread+ilen(lancuch)+1
2030 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2033 c----------------------------------------------------------------------------
2034 subroutine openunits
2036 include 'DIMENSIONS'
2040 character*16 form,nodename
2043 include 'COMMON.SETUP'
2044 include 'COMMON.IOUNITS'
2046 include 'COMMON.CONTROL'
2047 integer lenpre,lenpot,ilen,lentmp,npos
2049 character*3 out1file_text,ucase
2054 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2055 call getenv_loc("PREFIX",prefix)
2057 call getenv_loc("POT",pot)
2058 call getenv_loc("DIRTMP",tmpdir)
2059 call getenv_loc("CURDIR",curdir)
2060 call getenv_loc("OUT1FILE",out1file_text)
2061 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2062 out1file_text=ucase(out1file_text)
2063 if (out1file_text(1:1).eq."Y") then
2066 out1file=fg_rank.gt.0
2071 if (lentmp.gt.0) then
2072 write (*,'(80(1h!))')
2073 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2074 write (*,'(80(1h!))')
2075 write (*,*)"All output files will be on node /tmp directory."
2077 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2078 if (me.eq.king) then
2079 write (*,*) "The master node is ",nodename
2080 else if (fg_rank.eq.0) then
2081 write (*,*) "I am the CG slave node ",nodename
2083 write (*,*) "I am the FG slave node ",nodename
2086 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2087 lenpre = lentmp+lenpre+1
2089 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2090 C Get the names and open the input files
2091 #if defined(WINIFL) || defined(WINPGI)
2092 open(1,file=pref_orig(:ilen(pref_orig))//
2093 & '.inp',status='old',readonly,shared)
2094 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2095 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2096 C Get parameter filenames and open the parameter files.
2097 call getenv_loc('BONDPAR',bondname)
2098 open (ibond,file=bondname,status='old',readonly,shared)
2099 call getenv_loc('THETPAR',thetname)
2100 open (ithep,file=thetname,status='old',readonly,shared)
2101 call getenv_loc('ROTPAR',rotname)
2102 open (irotam,file=rotname,status='old',readonly,shared)
2103 call getenv_loc('TORPAR',torname)
2104 open (itorp,file=torname,status='old',readonly,shared)
2105 call getenv_loc('TORDPAR',tordname)
2106 open (itordp,file=tordname,status='old',readonly,shared)
2107 call getenv_loc('FOURIER',fouriername)
2108 open (ifourier,file=fouriername,status='old',readonly,shared)
2109 call getenv_loc('ELEPAR',elename)
2110 open (ielep,file=elename,status='old',readonly,shared)
2111 call getenv_loc('SIDEPAR',sidename)
2112 open (isidep,file=sidename,status='old',readonly,shared)
2113 call getenv_loc('LIPTRANPAR',liptranname)
2114 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2115 #elif (defined CRAY) || (defined AIX)
2116 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2118 c print *,"Processor",myrank," opened file 1"
2119 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2120 c print *,"Processor",myrank," opened file 9"
2121 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2122 C Get parameter filenames and open the parameter files.
2123 call getenv_loc('BONDPAR',bondname)
2124 open (ibond,file=bondname,status='old',action='read')
2125 c print *,"Processor",myrank," opened file IBOND"
2126 call getenv_loc('THETPAR',thetname)
2127 open (ithep,file=thetname,status='old',action='read')
2129 call getenv_loc('THETPARPDB',thetname_pdb)
2130 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2132 c print *,"Processor",myrank," opened file ITHEP"
2133 call getenv_loc('ROTPAR',rotname)
2134 open (irotam,file=rotname,status='old',action='read')
2136 call getenv_loc('ROTPARPDB',rotname_pdb)
2137 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2139 c print *,"Processor",myrank," opened file IROTAM"
2140 call getenv_loc('TORPAR',torname)
2141 open (itorp,file=torname,status='old',action='read')
2142 c print *,"Processor",myrank," opened file ITORP"
2143 call getenv_loc('TORDPAR',tordname)
2144 open (itordp,file=tordname,status='old',action='read')
2145 c print *,"Processor",myrank," opened file ITORDP"
2146 call getenv_loc('SCCORPAR',sccorname)
2147 open (isccor,file=sccorname,status='old',action='read')
2148 c print *,"Processor",myrank," opened file ISCCOR"
2149 call getenv_loc('FOURIER',fouriername)
2150 open (ifourier,file=fouriername,status='old',action='read')
2151 c print *,"Processor",myrank," opened file IFOURIER"
2152 call getenv_loc('ELEPAR',elename)
2153 open (ielep,file=elename,status='old',action='read')
2154 c print *,"Processor",myrank," opened file IELEP"
2155 call getenv_loc('SIDEPAR',sidename)
2156 open (isidep,file=sidename,status='old',action='read')
2157 call getenv_loc('LIPTRANPAR',liptranname)
2158 open (iliptranpar,file=liptranname,status='old',action='read')
2159 c print *,"Processor",myrank," opened file ISIDEP"
2160 c print *,"Processor",myrank," opened parameter files"
2162 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2163 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2164 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2165 C Get parameter filenames and open the parameter files.
2166 call getenv_loc('BONDPAR',bondname)
2167 open (ibond,file=bondname,status='old')
2168 call getenv_loc('THETPAR',thetname)
2169 open (ithep,file=thetname,status='old')
2171 call getenv_loc('THETPARPDB',thetname_pdb)
2172 open (ithep_pdb,file=thetname_pdb,status='old')
2174 call getenv_loc('ROTPAR',rotname)
2175 open (irotam,file=rotname,status='old')
2177 call getenv_loc('ROTPARPDB',rotname_pdb)
2178 open (irotam_pdb,file=rotname_pdb,status='old')
2180 call getenv_loc('TORPAR',torname)
2181 open (itorp,file=torname,status='old')
2182 call getenv_loc('TORDPAR',tordname)
2183 open (itordp,file=tordname,status='old')
2184 call getenv_loc('SCCORPAR',sccorname)
2185 open (isccor,file=sccorname,status='old')
2186 call getenv_loc('FOURIER',fouriername)
2187 open (ifourier,file=fouriername,status='old')
2188 call getenv_loc('ELEPAR',elename)
2189 open (ielep,file=elename,status='old')
2190 call getenv_loc('SIDEPAR',sidename)
2191 open (isidep,file=sidename,status='old')
2192 call getenv_loc('LIPTRANPAR',liptranname)
2193 open (iliptranpar,file=liptranname,status='old')
2195 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2197 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2198 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2199 C Get parameter filenames and open the parameter files.
2200 call getenv_loc('BONDPAR',bondname)
2201 open (ibond,file=bondname,status='old',readonly)
2202 call getenv_loc('THETPAR',thetname)
2203 open (ithep,file=thetname,status='old',readonly)
2204 call getenv_loc('ROTPAR',rotname)
2205 open (irotam,file=rotname,status='old',readonly)
2206 call getenv_loc('TORPAR',torname)
2207 open (itorp,file=torname,status='old',readonly)
2208 call getenv_loc('TORDPAR',tordname)
2209 open (itordp,file=tordname,status='old',readonly)
2210 call getenv_loc('SCCORPAR',sccorname)
2211 open (isccor,file=sccorname,status='old',readonly)
2213 call getenv_loc('THETPARPDB',thetname_pdb)
2214 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2216 call getenv_loc('FOURIER',fouriername)
2217 open (ifourier,file=fouriername,status='old',readonly)
2218 call getenv_loc('ELEPAR',elename)
2219 open (ielep,file=elename,status='old',readonly)
2220 call getenv_loc('SIDEPAR',sidename)
2221 open (isidep,file=sidename,status='old',readonly)
2222 call getenv_loc('LIPTRANPAR',liptranname)
2223 open (iliptranpar,file=liptranname,status='old',action='read')
2225 call getenv_loc('ROTPARPDB',rotname_pdb)
2226 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2230 call getenv_loc('TUBEPAR',tubename)
2231 #if defined(WINIFL) || defined(WINPGI)
2232 open (itube,file=tubename,status='old',readonly,shared)
2233 #elif (defined CRAY) || (defined AIX)
2234 open (itube,file=tubename,status='old',action='read')
2236 open (itube,file=tubename,status='old')
2238 open (itube,file=tubename,status='old',readonly)
2243 C 8/9/01 In the newest version SCp interaction constants are read from a file
2244 C Use -DOLDSCP to use hard-coded constants instead.
2246 call getenv_loc('SCPPAR',scpname)
2247 #if defined(WINIFL) || defined(WINPGI)
2248 open (iscpp,file=scpname,status='old',readonly,shared)
2249 #elif (defined CRAY) || (defined AIX)
2250 open (iscpp,file=scpname,status='old',action='read')
2252 open (iscpp,file=scpname,status='old')
2254 open (iscpp,file=scpname,status='old',readonly)
2257 call getenv_loc('PATTERN',patname)
2258 #if defined(WINIFL) || defined(WINPGI)
2259 open (icbase,file=patname,status='old',readonly,shared)
2260 #elif (defined CRAY) || (defined AIX)
2261 open (icbase,file=patname,status='old',action='read')
2263 open (icbase,file=patname,status='old')
2265 open (icbase,file=patname,status='old',readonly)
2268 C Open output file only for CG processes
2269 c print *,"Processor",myrank," fg_rank",fg_rank
2270 if (fg_rank.eq.0) then
2272 if (nodes.eq.1) then
2275 npos = dlog10(dfloat(nodes-1))+1
2277 if (npos.lt.3) npos=3
2278 write (liczba,'(i1)') npos
2279 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2281 write (liczba,form) me
2282 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2283 & liczba(:ilen(liczba))
2284 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2286 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2288 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2289 & liczba(:ilen(liczba))//'.mol2'
2290 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2291 & liczba(:ilen(liczba))//'.stat'
2293 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2294 & //liczba(:ilen(liczba))//'.stat')
2295 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2298 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2299 & liczba(:ilen(liczba))//'.const'
2304 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2305 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2306 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2307 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2308 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2310 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2312 rest2name=prefix(:ilen(prefix))//'.rst'
2314 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2317 #if defined(AIX) || defined(PGI) || defined(CRAY)
2318 if (me.eq.king .or. .not. out1file) then
2319 open(iout,file=outname,status='unknown')
2321 open(iout,file="/dev/null",status="unknown")
2325 if (fg_rank.gt.0) then
2326 write (liczba,'(i3.3)') myrank/nfgtasks
2327 write (ll,'(bz,i3.3)') fg_rank
2328 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2334 open(igeom,file=intname,status='unknown',position='append')
2335 open(ipdb,file=pdbname,status='unknown')
2336 open(imol2,file=mol2name,status='unknown')
2337 open(istat,file=statname,status='unknown',position='append')
2339 c1out open(iout,file=outname,status='unknown')
2342 if (me.eq.king .or. .not.out1file)
2343 & open(iout,file=outname,status='unknown')
2345 if (fg_rank.gt.0) then
2346 write (liczba,'(i3.3)') myrank/nfgtasks
2347 write (ll,'(bz,i3.3)') fg_rank
2348 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2353 open(igeom,file=intname,status='unknown',access='append')
2354 open(ipdb,file=pdbname,status='unknown')
2355 open(imol2,file=mol2name,status='unknown')
2356 open(istat,file=statname,status='unknown',access='append')
2358 c1out open(iout,file=outname,status='unknown')
2361 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2362 csa_seed=prefix(:lenpre)//'.CSA.seed'
2363 csa_history=prefix(:lenpre)//'.CSA.history'
2364 csa_bank=prefix(:lenpre)//'.CSA.bank'
2365 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2366 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2367 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2368 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2369 csa_int=prefix(:lenpre)//'.int'
2370 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2371 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2372 csa_in=prefix(:lenpre)//'.CSA.in'
2373 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2376 write (iout,'(80(1h-))')
2377 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2378 write (iout,'(80(1h-))')
2379 write (iout,*) "Input file : ",
2380 & pref_orig(:ilen(pref_orig))//'.inp'
2381 write (iout,*) "Output file : ",
2382 & outname(:ilen(outname))
2384 write (iout,*) "Sidechain potential file : ",
2385 & sidename(:ilen(sidename))
2387 write (iout,*) "SCp potential file : ",
2388 & scpname(:ilen(scpname))
2390 write (iout,*) "Electrostatic potential file : ",
2391 & elename(:ilen(elename))
2392 write (iout,*) "Cumulant coefficient file : ",
2393 & fouriername(:ilen(fouriername))
2394 write (iout,*) "Torsional parameter file : ",
2395 & torname(:ilen(torname))
2396 write (iout,*) "Double torsional parameter file : ",
2397 & tordname(:ilen(tordname))
2398 write (iout,*) "SCCOR parameter file : ",
2399 & sccorname(:ilen(sccorname))
2400 write (iout,*) "Bond & inertia constant file : ",
2401 & bondname(:ilen(bondname))
2402 write (iout,*) "Bending-torsion parameter file : ",
2403 & thetname(:ilen(thetname))
2404 write (iout,*) "Rotamer parameter file : ",
2405 & rotname(:ilen(rotname))
2406 write (iout,*) "Thetpdb parameter file : ",
2407 & thetname_pdb(:ilen(thetname_pdb))
2408 write (iout,*) "Threading database : ",
2409 & patname(:ilen(patname))
2411 &write (iout,*)" DIRTMP : ",
2413 write (iout,'(80(1h-))')
2417 c----------------------------------------------------------------------------
2418 subroutine card_concat(card)
2419 implicit real*8 (a-h,o-z)
2420 include 'DIMENSIONS'
2421 include 'COMMON.IOUNITS'
2423 character*80 karta,ucase
2425 read (inp,'(a)') karta
2428 do while (karta(80:80).eq.'&')
2429 card=card(:ilen(card)+1)//karta(:79)
2430 read (inp,'(a)') karta
2433 card=card(:ilen(card)+1)//karta
2436 c------------------------------------------------------------------------------
2439 include 'DIMENSIONS'
2440 include 'COMMON.CHAIN'
2441 include 'COMMON.IOUNITS'
2442 include 'COMMON.CONTROL'
2444 include 'COMMON.QRESTR'
2446 include 'COMMON.LAGRANGE.5diag'
2448 include 'COMMON.LAGRANGE'
2451 open(irest2,file=rest2name,status='unknown')
2452 read(irest2,*) totT,EK,potE,totE,t_bath
2455 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2458 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2461 read (irest2,*) iset
2466 c------------------------------------------------------------------------------
2467 subroutine read_fragments
2469 include 'DIMENSIONS'
2473 include 'COMMON.SETUP'
2474 include 'COMMON.CHAIN'
2475 include 'COMMON.IOUNITS'
2477 include 'COMMON.QRESTR'
2478 include 'COMMON.CONTROL'
2480 read(inp,*) nset,nfrag,npair,nfrag_back
2481 loc_qlike=(nfrag_back.lt.0)
2482 nfrag_back=iabs(nfrag_back)
2483 if(me.eq.king.or..not.out1file)
2484 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2485 & " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2487 read(inp,*) mset(iset)
2489 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2491 if(me.eq.king.or..not.out1file)
2492 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2493 & ifrag(2,i,iset), qinfrag(i,iset)
2496 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2498 if(me.eq.king.or..not.out1file)
2499 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2500 & ipair(2,i,iset), qinpair(i,iset)
2504 read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2505 & wfrag_back(2,i,iset),qin_back(2,i,iset),
2506 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2507 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2508 if(me.eq.king.or..not.out1file)
2509 & write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2510 & wfrag_back(2,i,iset),qin_back(3,i,iset),
2511 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2512 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2516 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2517 & wfrag_back(3,i,iset),
2518 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2519 if(me.eq.king.or..not.out1file)
2520 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2521 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2527 C---------------------------------------------------------------------------
2528 subroutine read_afminp
2530 include 'DIMENSIONS'
2534 include 'COMMON.SETUP'
2535 include 'COMMON.CONTROL'
2536 include 'COMMON.CHAIN'
2537 include 'COMMON.IOUNITS'
2538 include 'COMMON.SBRIDGE'
2539 character*320 afmcard
2541 c print *, "wchodze"
2542 call card_concat(afmcard)
2543 call readi(afmcard,"BEG",afmbeg,0)
2544 call readi(afmcard,"END",afmend,0)
2545 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2546 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2547 c print *,'FORCE=' ,forceAFMconst
2548 CCCC NOW PROPERTIES FOR AFM
2551 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2553 distafminit=dsqrt(distafminit)
2554 c print *,'initdist',distafminit
2557 c-------------------------------------------------------------------------------
2558 subroutine read_saxs_constr
2560 include 'DIMENSIONS'
2564 include 'COMMON.SETUP'
2565 include 'COMMON.CONTROL'
2566 include 'COMMON.SAXS'
2567 include 'COMMON.CHAIN'
2568 include 'COMMON.IOUNITS'
2569 include 'COMMON.SBRIDGE'
2570 double precision cm(3),cnorm
2573 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2575 if (saxs_mode.eq.0) then
2576 c SAXS distance distribution
2578 read(inp,*) distsaxs(i),Psaxs(i)
2582 Cnorm = Cnorm + Psaxs(i)
2584 write (iout,*) "Cnorm",Cnorm
2586 Psaxs(i)=Psaxs(i)/Cnorm
2588 write (iout,*) "Normalized distance distribution from SAXS"
2590 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2594 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2596 write (iout,*) "Wsaxs0",Wsaxs0
2600 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2607 cm(j)=cm(j)+Csaxs(j,i)
2615 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2618 write (iout,*) "SAXS sphere coordinates"
2620 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2626 c-------------------------------------------------------------------------------
2627 subroutine read_dist_constr
2629 include 'DIMENSIONS'
2633 include 'COMMON.SETUP'
2634 include 'COMMON.CONTROL'
2635 include 'COMMON.CHAIN'
2636 include 'COMMON.IOUNITS'
2637 include 'COMMON.SBRIDGE'
2638 include 'COMMON.INTERACT'
2639 integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
2640 integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
2641 double precision wfrag_(100),wpair_(1000)
2642 double precision ddjk,dist,dist_cut,fordepthmax
2643 character*5000 controlcard
2644 logical normalize,next
2646 double precision scal_bfac
2647 double precision xlink(4,0:4) /
2649 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2650 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2651 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2652 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2653 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2654 c print *, "WCHODZE"
2655 write (iout,*) "Calling read_dist_constr"
2656 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2658 restr_on_coord=.false.
2667 call card_concat(controlcard)
2668 next = index(controlcard,"NEXT").gt.0
2669 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2670 write (iout,*) "restr_type",restr_type
2671 call readi(controlcard,"NFRAG",nfrag_,0)
2672 call readi(controlcard,"NFRAG",nfrag_,0)
2673 call readi(controlcard,"NPAIR",npair_,0)
2674 call readi(controlcard,"NDIST",ndist_,0)
2675 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2676 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2677 if (restr_type.eq.10)
2678 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2679 if (restr_type.eq.12)
2680 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2681 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2682 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2683 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2684 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2685 normalize = index(controlcard,"NORMALIZE").gt.0
2686 write (iout,*) "WBOLTZD",wboltzd
2687 write (iout,*) "SCAL_PEAK",scal_peak
2688 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2689 write (iout,*) "IFRAG"
2691 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2693 write (iout,*) "IPAIR"
2695 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2697 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2699 & "Distance restraints as generated from reference structure"
2701 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2702 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2703 & ifrag_(2,i)=nstart_sup+nsup-1
2704 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2706 if (wfrag_(i).eq.0.0d0) cycle
2707 do j=ifrag_(1,i),ifrag_(2,i)-1
2708 do k=j+1,ifrag_(2,i)
2709 c write (iout,*) "j",j," k",k
2711 if (restr_type.eq.1) then
2717 forcon(nhpb)=wfrag_(i)
2718 else if (constr_dist.eq.2) then
2719 if (ddjk.le.dist_cut) then
2725 forcon(nhpb)=wfrag_(i)
2727 else if (restr_type.eq.3) then
2733 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2736 if (.not.out1file .or. me.eq.king)
2737 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2738 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2740 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2741 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2747 if (wpair_(i).eq.0.0d0) cycle
2755 do j=ifrag_(1,ii),ifrag_(2,ii)
2756 do k=ifrag_(1,jj),ifrag_(2,jj)
2758 if (restr_type.eq.1) then
2764 forcon(nhpb)=wpair_(i)
2765 else if (constr_dist.eq.2) then
2766 if (ddjk.le.dist_cut) then
2772 forcon(nhpb)=wpair_(i)
2774 else if (restr_type.eq.3) then
2780 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2783 if (.not.out1file .or. me.eq.king)
2784 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2785 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2787 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2788 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2795 write (iout,*) "Distance restraints as read from input"
2797 if (restr_type.eq.12) then
2798 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2799 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2800 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2801 & fordepth_peak(nhpb_peak+1),npeak
2802 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2803 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2804 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2805 c & fordepth_peak(nhpb_peak+1),npeak
2806 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2807 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2808 nhpb_peak=nhpb_peak+1
2809 irestr_type_peak(nhpb_peak)=12
2810 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2813 if (.not.out1file .or. me.eq.king)
2814 & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2815 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2816 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2817 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2818 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2820 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2821 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2822 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2823 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2824 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2826 if (ibecarb_peak(nhpb_peak).eq.3) then
2827 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2828 else if (ibecarb_peak(nhpb_peak).eq.2) then
2829 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2830 else if (ibecarb_peak(nhpb_peak).eq.1) then
2831 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2832 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2834 else if (restr_type.eq.11) then
2835 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2836 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2837 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2838 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2840 irestr_type(nhpb)=11
2842 if (.not.out1file .or. me.eq.king)
2843 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2844 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2845 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2847 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2848 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2849 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2851 c if (ibecarb(nhpb).gt.0) then
2852 c ihpb(nhpb)=ihpb(nhpb)+nres
2853 c jhpb(nhpb)=jhpb(nhpb)+nres
2855 if (ibecarb(nhpb).eq.3) then
2856 ihpb(nhpb)=ihpb(nhpb)+nres
2857 else if (ibecarb(nhpb).eq.2) then
2858 ihpb(nhpb)=ihpb(nhpb)+nres
2859 else if (ibecarb(nhpb).eq.1) then
2860 ihpb(nhpb)=ihpb(nhpb)+nres
2861 jhpb(nhpb)=jhpb(nhpb)+nres
2863 else if (restr_type.eq.10) then
2864 c Cross-lonk Markov-like potential
2865 call card_concat(controlcard)
2866 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2867 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2869 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2870 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2871 if (index(controlcard,"ZL").gt.0) then
2873 else if (index(controlcard,"ADH").gt.0) then
2875 else if (index(controlcard,"PDH").gt.0) then
2877 else if (index(controlcard,"DSS").gt.0) then
2882 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2883 & xlink(1,link_type))
2884 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2885 & xlink(2,link_type))
2886 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2887 & xlink(3,link_type))
2888 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2889 & xlink(4,link_type))
2890 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2891 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2892 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2893 if (forcon(nhpb+1).le.0.0d0 .or.
2894 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2896 irestr_type(nhpb)=10
2897 if (ibecarb(nhpb).eq.3) then
2898 jhpb(nhpb)=jhpb(nhpb)+nres
2899 else if (ibecarb(nhpb).eq.2) then
2900 ihpb(nhpb)=ihpb(nhpb)+nres
2901 else if (ibecarb(nhpb).eq.1) then
2902 ihpb(nhpb)=ihpb(nhpb)+nres
2903 jhpb(nhpb)=jhpb(nhpb)+nres
2906 if (.not.out1file .or. me.eq.king)
2907 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2908 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2909 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2912 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2913 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2914 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2919 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2920 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2921 if (forcon(nhpb+1).gt.0.0d0) then
2923 if (dhpb1(nhpb).eq.0.0d0) then
2928 if (ibecarb(nhpb).eq.3) then
2929 jhpb(nhpb)=jhpb(nhpb)+nres
2930 else if (ibecarb(nhpb).eq.2) then
2931 ihpb(nhpb)=ihpb(nhpb)+nres
2932 else if (ibecarb(nhpb).eq.1) then
2933 ihpb(nhpb)=ihpb(nhpb)+nres
2934 jhpb(nhpb)=jhpb(nhpb)+nres
2936 if (dhpb(nhpb).eq.0.0d0)
2937 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2940 if (.not.out1file .or. me.eq.king)
2941 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2942 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2944 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2945 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2948 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2949 C if (forcon(nhpb+1).gt.0.0d0) then
2951 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2954 if (restr_type.eq.4) then
2955 write (iout,*) "The BFAC array"
2957 write (iout,'(i5,f10.5)') i,bfac(i)
2960 if (itype(i).eq.ntyp1) cycle
2962 if (itype(j).eq.ntyp1) cycle
2963 if (itype(i).eq.10) then
2968 if (itype(j).eq.10) then
2978 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2981 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2982 ihpb(nhpb)=i+nres*ii
2983 jhpb(nhpb)=j+nres*jj
2984 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2986 if (.not.out1file .or. me.eq.king) then
2987 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2988 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2989 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2993 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2994 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2995 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
3005 if (restr_type.eq.5) then
3006 restr_on_coord=.true.
3008 if (itype(i).eq.ntyp1) cycle
3009 bfac(i)=(scal_bfac/bfac(i))**2
3018 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
3019 & fordepthmax=fordepth(i)
3022 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
3027 c-------------------------------------------------------------------------------
3028 subroutine read_constr_homology
3030 include 'DIMENSIONS'
3034 include 'COMMON.SETUP'
3035 include 'COMMON.CONTROL'
3036 include 'COMMON.HOMOLOGY'
3037 include 'COMMON.CHAIN'
3038 include 'COMMON.IOUNITS'
3040 include 'COMMON.QRESTR'
3041 include 'COMMON.GEO'
3042 include 'COMMON.INTERACT'
3043 include 'COMMON.NAMES'
3045 c For new homol impl
3047 include 'COMMON.VAR'
3050 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
3052 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
3053 c & sigma_odl_temp(maxres,maxres,max_template)
3055 character*24 model_ki_dist, model_ki_angle
3056 character*500 controlcard
3057 integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
3058 & ik,iistart,nres_temp
3061 logical liiflag,lfirst
3064 c FP - Nov. 2014 Temporary specifications for new vars
3066 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
3068 double precision, dimension (max_template,maxres) :: rescore
3069 double precision, dimension (max_template,maxres) :: rescore2
3070 double precision, dimension (max_template,maxres) :: rescore3
3071 double precision distal
3072 character*24 tpl_k_rescore
3073 character*256 pdbfile
3074 c -----------------------------------------------------------------
3075 c Reading multiple PDB ref structures and calculation of retraints
3076 c not using pre-computed ones stored in files model_ki_{dist,angle}
3078 c -----------------------------------------------------------------
3081 c Alternative: reading from input
3082 call card_concat(controlcard)
3083 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
3084 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
3085 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
3086 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
3087 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3088 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3089 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3090 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3091 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3092 if(.not.read2sigma.and.start_from_model) then
3093 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3094 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3095 start_from_model=.false.
3096 iranconf=(indpdb.le.0)
3098 if(start_from_model .and. (me.eq.king .or. .not. out1file))
3099 & write(iout,*) 'START_FROM_MODELS is ON'
3100 c if(start_from_model .and. rest) then
3101 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3102 c write(iout,*) 'START_FROM_MODELS is OFF'
3103 c write(iout,*) 'remove restart keyword from input'
3106 if (start_from_model) nmodel_start=constr_homology
3107 if (homol_nset.gt.1)then
3108 call card_concat(controlcard)
3109 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3110 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3111 write(iout,*) "iset homology_weight "
3113 write(iout,*) i,waga_homology(i)
3116 iset=mod(kolor,homol_nset)+1
3119 waga_homology(1)=1.0
3122 cd write (iout,*) "nnt",nnt," nct",nct
3129 c write(iout,*) 'nnt=',nnt,'nct=',nct
3132 do k=1,constr_homology
3145 if (read_homol_frag) then
3146 call read_klapaucjusz
3149 do k=1,constr_homology
3151 read(inp,'(a)') pdbfile
3152 pdbfiles_chomo(k)=pdbfile
3153 if(me.eq.king .or. .not. out1file)
3154 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3155 & pdbfile(:ilen(pdbfile))
3156 open(ipdbin,file=pdbfile,status='old',err=33)
3158 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3159 & pdbfile(:ilen(pdbfile))
3162 c print *,'Begin reading pdb data'
3164 c Files containing res sim or local scores (former containing sigmas)
3167 write(kic2,'(bz,i2.2)') k
3169 tpl_k_rescore="template"//kic2//".sco"
3173 if (read2sigma) then
3174 call readpdb_template(k)
3181 c Distance restraints
3184 C Copy the coordinates from reference coordinates (?)
3185 do i=1,2*nres_chomo(k)
3188 c write (iout,*) "c(",j,i,") =",c(j,i)
3192 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3194 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3195 open (ientin,file=tpl_k_rescore,status='old')
3196 if (nnt.gt.1) rescore(k,1)=0.0d0
3197 do irec=nnt,nct ! loop for reading res sim
3198 if (read2sigma) then
3199 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3200 & rescore3_tmp,idomain_tmp
3202 idomain(k,i_tmp)=idomain_tmp
3203 rescore(k,i_tmp)=rescore_tmp
3204 rescore2(k,i_tmp)=rescore2_tmp
3205 rescore3(k,i_tmp)=rescore3_tmp
3206 if (.not. out1file .or. me.eq.king)
3207 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3208 & i_tmp,rescore2_tmp,rescore_tmp,
3209 & rescore3_tmp,idomain_tmp
3212 read (ientin,*,end=1401) rescore_tmp
3214 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3215 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3216 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3221 if (waga_dist.ne.0.0d0) then
3229 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3230 c write (iout,*) k,i,j,distal,dist2_cut
3232 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3233 & .and. distal.le.dist2_cut ) then
3239 c write (iout,*) "k",k
3240 c write (iout,*) "i",i," j",j," constr_homology",
3245 if (read2sigma) then
3248 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3250 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3251 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3252 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3254 if (odl(k,ii).le.dist_cut) then
3255 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3258 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3259 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3261 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3262 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3266 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3269 l_homo(k,ii)=.false.
3275 c write (iout,*) "Distance restraints set"
3278 c Theta, dihedral and SC retraints
3280 if (waga_angle.gt.0.0d0) then
3281 c open (ientin,file=tpl_k_sigma_dih,status='old')
3282 c do irec=1,maxres-3 ! loop for reading sigma_dih
3283 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3284 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3285 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3286 c & sigma_dih(k,i+nnt-1)
3291 if (idomain(k,i).eq.0) then
3295 dih(k,i)=phiref(i) ! right?
3296 c read (ientin,*) sigma_dih(k,i) ! original variant
3297 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3298 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3299 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3300 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3301 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3303 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3304 & rescore(k,i-2)+rescore(k,i-3))/4.0
3305 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3306 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3307 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3308 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3309 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3310 c Instead of res sim other local measure of b/b str reliability possible
3311 if (sigma_dih(k,i).ne.0)
3312 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3313 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3317 c write (iout,*) "Dihedral angle restraints set"
3320 if (waga_theta.gt.0.0d0) then
3321 c open (ientin,file=tpl_k_sigma_theta,status='old')
3322 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3323 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3324 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3325 c & sigma_theta(k,i+nnt-1)
3330 do i = nnt+2,nct ! right? without parallel.
3331 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3332 c do i=ithet_start,ithet_end ! with FG parallel.
3333 if (idomain(k,i).eq.0) then
3334 sigma_theta(k,i)=0.0
3337 thetatpl(k,i)=thetaref(i)
3338 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3339 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3340 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3341 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3342 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3343 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3344 & rescore(k,i-2))/3.0
3345 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3346 if (sigma_theta(k,i).ne.0)
3347 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3349 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3350 c rescore(k,i-2) ! right expression ?
3351 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3354 c write (iout,*) "Angle restraints set"
3357 if (waga_d.gt.0.0d0) then
3358 c open (ientin,file=tpl_k_sigma_d,status='old')
3359 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3360 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3361 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3362 c & sigma_d(k,i+nnt-1)
3366 do i = nnt,nct ! right? without parallel.
3367 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3368 c do i=loc_start,loc_end ! with FG parallel.
3369 if (itype(i).eq.10) cycle
3370 if (idomain(k,i).eq.0 ) then
3377 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3378 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3379 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3380 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3381 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3382 if (sigma_d(k,i).ne.0)
3383 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3385 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3386 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3387 c read (ientin,*) sigma_d(k,i) ! 1st variant
3391 c write (iout,*) "SC restraints set"
3394 c remove distance restraints not used in any model from the list
3395 c shift data in all arrays
3397 c write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
3398 if (waga_dist.ne.0.0d0) then
3405 c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3406 c & .and. distal.le.dist2_cut ) then
3407 c write (iout,*) "i",i," j",j," ii",ii
3409 if (ii_in_use(ii).eq.0.and.liiflag.or.
3410 & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3417 if(i10.eq.lim_odl) i10=i10+1
3419 ires_homo(iistart+ki)=ires_homo(ki+i01)
3420 jres_homo(iistart+ki)=jres_homo(ki+i01)
3421 ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3422 do k=1,constr_homology
3423 odl(k,iistart+ki)=odl(k,ki+i01)
3424 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3425 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3428 iistart=iistart+i10-i01
3431 if (ii_in_use(ii).ne.0.and..not.liiflag) then
3439 c write (iout,*) "Removing distances completed"
3441 endif ! .not. klapaucjusz
3443 if (constr_homology.gt.0) call homology_partition
3444 c write (iout,*) "After homology_partition"
3446 if (constr_homology.gt.0) call init_int_table
3447 c write (iout,*) "After init_int_table"
3449 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3450 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3454 if (.not.out_template_restr) return
3455 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3456 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3457 write (iout,*) "Distance restraints from templates"
3459 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3460 & ii,ires_homo(ii),jres_homo(ii),
3461 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3462 & ki=1,constr_homology)
3464 write (iout,*) "Dihedral angle restraints from templates"
3466 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3467 & (rad2deg*dih(ki,i),
3468 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3470 write (iout,*) "Virtual-bond angle restraints from templates"
3472 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3473 & (rad2deg*thetatpl(ki,i),
3474 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3476 write (iout,*) "SC restraints from templates"
3478 write(iout,'(i5,100(4f8.2,4x))') i,
3479 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3480 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3483 c -----------------------------------------------------------------
3486 c----------------------------------------------------------------------
3488 subroutine flush(iu)
3493 subroutine flush(iu)
3498 c------------------------------------------------------------------------------
3499 subroutine copy_to_tmp(source)
3501 include "DIMENSIONS"
3502 include "COMMON.IOUNITS"
3503 character*(*) source
3504 character* 256 tmpfile
3508 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3509 inquire(file=tmpfile,exist=ex)
3511 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3512 & " to temporary directory..."
3513 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3514 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3518 c------------------------------------------------------------------------------
3519 subroutine move_from_tmp(source)
3521 include "DIMENSIONS"
3522 include "COMMON.IOUNITS"
3523 character*(*) source
3526 write (*,*) "Moving ",source(:ilen(source)),
3527 & " from temporary directory to working directory"
3528 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3529 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3532 c------------------------------------------------------------------------------
3533 subroutine random_init(seed)
3535 C Initialize random number generator
3538 include 'DIMENSIONS'
3541 logical OKRandom, prng_restart
3543 integer iseed_array(4)
3544 integer error_msg,ierr
3546 include 'COMMON.IOUNITS'
3547 include 'COMMON.TIME1'
3548 include 'COMMON.THREAD'
3549 include 'COMMON.SBRIDGE'
3550 include 'COMMON.CONTROL'
3551 include 'COMMON.MCM'
3552 include 'COMMON.MAP'
3553 include 'COMMON.HEADER'
3554 include 'COMMON.CSA'
3555 include 'COMMON.CHAIN'
3556 include 'COMMON.MUCA'
3558 include 'COMMON.FFIELD'
3559 include 'COMMON.SETUP'
3561 double precision seed,ran_number
3562 iseed=-dint(dabs(seed))
3563 if (iseed.eq.0) then
3564 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3565 & 'Random seed undefined. The program will stop.'
3566 write (*,'(/80(1h*)/20x,a/80(1h*))')
3567 & 'Random seed undefined. The program will stop.'
3569 call mpi_finalize(mpi_comm_world,ierr)
3571 stop 'Bad random seed.'
3574 if (fg_rank.eq.0) then
3578 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3579 OKRandom = prng_restart(me,iseed)
3582 tmp=65536.0d0**(4-i)
3583 iseed_array(i) = dint(seed/tmp)
3584 seed=seed-iseed_array(i)*tmp
3587 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3588 & (iseed_array(i),i=1,4)
3589 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3590 & (iseed_array(i),i=1,4)
3591 OKRandom = prng_restart(me,iseed_array)
3594 c r1 = prng_next(me)
3595 r1=ran_number(0.0D0,1.0D0)
3597 & write (iout,*) 'ran_num',r1
3598 if (r1.lt.0.0d0) OKRandom=.false.
3600 if (.not.OKRandom) then
3601 write (iout,*) 'PRNG IS NOT WORKING!!!'
3602 print *,'PRNG IS NOT WORKING!!!'
3605 call mpi_abort(mpi_comm_world,error_msg,ierr)
3608 write (iout,*) 'too many processors for parallel prng'
3609 write (*,*) 'too many processors for parallel prng'
3617 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3621 c----------------------------------------------------------------------
3622 subroutine read_klapaucjusz
3624 include 'DIMENSIONS'
3628 include 'COMMON.SETUP'
3629 include 'COMMON.CONTROL'
3630 include 'COMMON.HOMOLOGY'
3631 include 'COMMON.CHAIN'
3632 include 'COMMON.IOUNITS'
3634 include 'COMMON.GEO'
3635 include 'COMMON.INTERACT'
3636 include 'COMMON.NAMES'
3637 character*256 fragfile
3638 integer ninclust(maxclust),inclust(max_template,maxclust),
3639 & nresclust(maxclust),iresclust(maxres,maxclust),nclust
3642 character*24 model_ki_dist, model_ki_angle
3643 character*500 controlcard
3644 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
3645 & ik,ll,ii,kk,iistart,iishift,lim_xx
3646 double precision distal
3647 logical lprn /.true./
3654 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3655 double precision, dimension (max_template,maxres) :: rescore
3656 double precision, dimension (max_template,maxres) :: rescore2
3657 character*24 tpl_k_rescore
3658 character*256 pdbfile
3661 c For new homol impl
3663 include 'COMMON.VAR'
3665 call getenv("FRAGFILE",fragfile)
3666 open(ientin,file=fragfile,status="old",err=10)
3667 read(ientin,*) constr_homology,nclust
3668 nmodel_start=constr_homology
3674 do k=1,constr_homology
3675 read(ientin,'(a)') pdbfile
3676 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3677 & pdbfile(:ilen(pdbfile))
3678 pdbfiles_chomo(k)=pdbfile
3679 open(ipdbin,file=pdbfile,status='old',err=33)
3681 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3682 & pdbfile(:ilen(pdbfile))
3687 call readpdb_template(k)
3697 read(ientin,*) ninclust(i),nresclust(i)
3698 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3699 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3702 c Loop over clusters
3705 do ll = 1,ninclust(l)
3713 idomain(k,iresclust(i,l)+1) = 1
3715 idomain(k,iresclust(i,l)) = 1
3719 c Distance restraints
3722 C Copy the coordinates from reference coordinates (?)
3728 c write (iout,*) "c(",j,i,") =",c(j,i)
3731 call int_from_cart(.true.,.false.)
3732 call sc_loc_geom(.false.)
3734 thetaref(i)=theta(i)
3738 if (waga_dist.ne.0.0d0) then
3746 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3747 c write (iout,*) k,i,j,distal,dist2_cut
3749 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3750 & .and. distal.le.dist2_cut ) then
3756 c write (iout,*) "k",k
3757 c write (iout,*) "i",i," j",j," constr_homology",
3762 if (read2sigma) then
3765 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3767 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3768 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3769 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3771 if (odl(k,ii).le.dist_cut) then
3772 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3775 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3776 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3778 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3779 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3783 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3786 c l_homo(k,ii)=.false.
3793 c Theta, dihedral and SC retraints
3795 if (waga_angle.gt.0.0d0) then
3797 if (idomain(k,i).eq.0) then
3798 c sigma_dih(k,i)=0.0
3802 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3803 & rescore(k,i-2)+rescore(k,i-3))/4.0
3804 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3805 c & " sigma_dihed",sigma_dih(k,i)
3806 if (sigma_dih(k,i).ne.0)
3807 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3812 if (waga_theta.gt.0.0d0) then
3814 if (idomain(k,i).eq.0) then
3815 c sigma_theta(k,i)=0.0
3818 thetatpl(k,i)=thetaref(i)
3819 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3820 & rescore(k,i-2))/3.0
3821 if (sigma_theta(k,i).ne.0)
3822 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3826 if (waga_d.gt.0.0d0) then
3828 if (itype(i).eq.10) cycle
3829 if (idomain(k,i).eq.0 ) then
3836 sigma_d(k,i)=rescore(k,i)
3837 if (sigma_d(k,i).ne.0)
3838 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3839 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3845 c remove distance restraints not used in any model from the list
3846 c shift data in all arrays
3848 if (waga_dist.ne.0.0d0) then
3854 if (ii_in_use(ii).eq.0.and.liiflag) then
3858 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3859 & .not.liiflag.and.ii.eq.lim_odl) then
3860 if (ii.eq.lim_odl) then
3861 iishift=ii-iistart+1
3866 do ki=iistart,lim_odl-iishift
3867 ires_homo(ki)=ires_homo(ki+iishift)
3868 jres_homo(ki)=jres_homo(ki+iishift)
3869 ii_in_use(ki)=ii_in_use(ki+iishift)
3870 do k=1,constr_homology
3871 odl(k,ki)=odl(k,ki+iishift)
3872 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3873 l_homo(k,ki)=l_homo(k,ki+iishift)
3877 lim_odl=lim_odl-iishift
3884 10 stop "Error in fragment file"