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 readi(controlcard,'PERMUT',npermut,1)
153 call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
154 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
155 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
156 c call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
157 c call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
158 c call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
159 c call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
160 c call reada(controlcard,'DRMS',drms,0.1D0)
161 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
162 c write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
163 c write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
164 c write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
165 c write (iout,'(a,f10.1)')'DRMS = ',drms
166 cc write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
167 c write (iout,'(a,f10.1)') 'Time limit (min):',timlim
169 call readi(controlcard,'NZ_START',nz_start,0)
170 call readi(controlcard,'NZ_END',nz_end,0)
171 call readi(controlcard,'NRAN_START',nran_start,0)
172 write (iout,*) "nran_start",nran_start
173 c call readi(controlcard,'IZ_SC',iz_sc,0)
175 safety = 60.0d0*safety
177 call readi(controlcard,"INTER_LIST_UPDATE",imatupdate,100)
178 call reada(controlcard,"T_BATH",t_bath,300.0d0)
179 minim=(index(controlcard,'MINIMIZE').gt.0)
180 dccart=(index(controlcard,'CART').gt.0)
181 overlapsc=(index(controlcard,'OVERLAP').gt.0)
182 overlapsc=.not.overlapsc
183 searchsc=(index(controlcard,'SEARCHSC').gt.0 .and.
184 & index(controlcard,'NOSEARCHSC').eq.0)
185 sideadd=(index(controlcard,'SIDEADD').gt.0)
186 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
187 mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
188 outpdb=(index(controlcard,'PDBOUT').gt.0)
189 outmol2=(index(controlcard,'MOL2OUT').gt.0)
190 pdbref=(index(controlcard,'PDBREF').gt.0)
191 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
192 indpdb=index(controlcard,'PDBSTART')
193 AFMlog=(index(controlcard,'AFM'))
194 selfguide=(index(controlcard,'SELFGUIDE'))
195 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
196 call readi(controlcard,'TUBEMOD',tubelog,0)
197 c write (iout,*) TUBElog,"TUBEMODE"
198 call readi(controlcard,'IPRINT',iprint,0)
199 C SHIELD keyword sets if the shielding effect of side-chains is used
200 C 0 denots no shielding is used all peptide are equally despite the
201 C solvent accesible area
202 C 1 the newly introduced function
203 C 2 reseved for further possible developement
204 call readi(controlcard,'SHIELD',shield_mode,0)
205 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
206 write(iout,*) "shield_mode",shield_mode
208 call readi(controlcard,'TORMODE',tor_mode,0)
209 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
210 write(iout,*) "torsional and valence angle mode",tor_mode
211 call readi(controlcard,'MAXGEN',maxgen,10000)
212 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
213 call readi(controlcard,"KDIAG",kdiag,0)
214 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
215 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
216 & write (iout,*) "RESCALE_MODE",rescale_mode
217 split_ene=index(controlcard,'SPLIT_ENE').gt.0
218 if (index(controlcard,'REGULAR').gt.0.0D0) then
219 call reada(controlcard,'WEIDIS',weidis,0.1D0)
224 if (index(controlcard,"CASC").gt.0) then
226 c else if (index(controlcard,"CAONLY").gt.0) then
228 else if (index(controlcard,"SCONLY").gt.0) then
232 c write (iout,*) "ERROR! Must specify CASC CAONLY or SCONLY when",
233 c & " specifying REFSTR, PDBREF or REGULAR."
237 if (index(controlcard,'CHECKGRAD').gt.0) then
239 if (index(controlcard,' CART').gt.0) then
241 elseif (index(controlcard,'CARINT').gt.0) then
246 call reada(controlcard,'DELTA',aincr,1.0d-7)
247 c write (iout,*) "icheckgrad",icheckgrad
248 elseif (index(controlcard,'THREAD').gt.0) then
250 call readi(controlcard,'THREAD',nthread,0)
251 if (nthread.gt.0) then
252 call reada(controlcard,'WEIDIS',weidis,0.1D0)
255 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
256 stop 'Error termination in Read_Control.'
258 else if (index(controlcard,'MCMA').gt.0) then
260 else if (index(controlcard,'MCEE').gt.0) then
262 else if (index(controlcard,'MULTCONF').gt.0) then
264 else if (index(controlcard,'MAP').gt.0) then
266 call readi(controlcard,'MAP',nmap,0)
267 else if (index(controlcard,'CSA').gt.0) then
269 crc else if (index(controlcard,'ZSCORE').gt.0) then
271 crc ZSCORE is rm from UNRES, modecalc=9 is available
274 cfcm else if (index(controlcard,'MCMF').gt.0) then
276 else if (index(controlcard,'SOFTREG').gt.0) then
278 else if (index(controlcard,'CHECK_BOND').gt.0) then
280 else if (index(controlcard,'TEST').gt.0) then
282 else if (index(controlcard,'MD').gt.0) then
284 else if (index(controlcard,'RE ').gt.0) then
288 lmuca=index(controlcard,'MUCA').gt.0
289 call readi(controlcard,'MUCADYN',mucadyn,0)
290 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
291 if (lmuca .and. (me.eq.king .or. .not.out1file ))
293 write (iout,*) 'MUCADYN=',mucadyn
294 write (iout,*) 'MUCASMOOTH=',muca_smooth
297 iscode=index(controlcard,'ONE_LETTER')
298 indphi=index(controlcard,'PHI')
299 indback=index(controlcard,'BACK')
300 iranconf=index(controlcard,'RAND_CONF')
301 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
302 extconf=(index(controlcard,'EXTCONF').gt.0)
303 if (start_from_model) then
307 i2ndstr=index(controlcard,'USE_SEC_PRED')
308 gradout=index(controlcard,'GRADOUT').gt.0
309 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
310 C DISTCHAINMAX become obsolete for periodic boundry condition
311 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
312 C Reading the dimensions of box in x,y,z coordinates
313 call reada(controlcard,'BOXX',boxxsize,100.0d0)
314 call reada(controlcard,'BOXY',boxysize,100.0d0)
315 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
316 write(iout,*) "Periodic box dimensions",boxxsize,boxysize,boxzsize
317 c Cutoff range for interactions
318 call reada(controlcard,"R_CUT_INT",r_cut_int,25.0d0)
319 call reada(controlcard,"R_CUT_RESPA",r_cut_respa,2.0d0)
320 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
321 write (iout,*) "Cutoff on interactions",r_cut_int
323 & "Cutoff in switching short and long range interactions in RESPA",
325 write (iout,*) "lambda in switch function",rlamb
326 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
327 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
328 if (lipthick.gt.0.0d0) then
329 bordliptop=(boxzsize+lipthick)/2.0
330 bordlipbot=bordliptop-lipthick
332 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
333 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
334 buflipbot=bordlipbot+lipbufthick
335 bufliptop=bordliptop-lipbufthick
336 if ((lipbufthick*2.0d0).gt.lipthick)
337 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
339 write(iout,*) "bordliptop=",bordliptop
340 write(iout,*) "bordlipbot=",bordlipbot
341 write(iout,*) "bufliptop=",bufliptop
342 write(iout,*) "buflipbot=",buflipbot
343 write (iout,*) "SHIELD MODE",shield_mode
344 if (TUBElog.gt.0) then
345 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
346 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
347 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
348 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
349 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
350 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
351 buftubebot=bordtubebot+tubebufthick
352 buftubetop=bordtubetop-tubebufthick
354 if (me.eq.king .or. .not.out1file )
355 & write (iout,*) "DISTCHAINMAX",distchainmax
357 if(me.eq.king.or..not.out1file)
358 & write (iout,'(2a)') diagmeth(kdiag),
359 & ' routine used to diagonalize matrices.'
362 c--------------------------------------------------------------------------
363 subroutine read_REMDpar
369 include 'COMMON.IOUNITS'
370 include 'COMMON.TIME1'
373 include 'COMMON.LANGEVIN'
376 include 'COMMON.LANGEVIN.lang0.5diag'
378 include 'COMMON.LANGEVIN.lang0'
381 include 'COMMON.INTERACT'
382 include 'COMMON.NAMES'
384 include 'COMMON.REMD'
385 include 'COMMON.CONTROL'
386 include 'COMMON.SETUP'
388 character*320 controlcard
389 character*3200 controlcard1
390 integer iremd_m_total,i
392 if(me.eq.king.or..not.out1file)
393 & write (iout,*) "REMD setup"
395 call card_concat(controlcard)
396 call readi(controlcard,"NREP",nrep,3)
397 call readi(controlcard,"NSTEX",nstex,1000)
398 call reada(controlcard,"RETMIN",retmin,10.0d0)
399 call reada(controlcard,"RETMAX",retmax,1000.0d0)
400 mremdsync=(index(controlcard,'SYNC').gt.0)
401 call readi(controlcard,"NSYN",i_sync_step,100)
402 restart1file=(index(controlcard,'REST1FILE').gt.0)
403 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
404 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
405 if(max_cache_traj_use.gt.max_cache_traj)
406 & max_cache_traj_use=max_cache_traj
407 if(me.eq.king.or..not.out1file) then
408 cd if (traj1file) then
409 crc caching is in testing - NTWX is not ignored
410 cd write (iout,*) "NTWX value is ignored"
411 cd write (iout,*) " trajectory is stored to one file by master"
412 cd write (iout,*) " before exchange at NSTEX intervals"
414 write (iout,*) "NREP= ",nrep
415 write (iout,*) "NSTEX= ",nstex
416 write (iout,*) "SYNC= ",mremdsync
417 write (iout,*) "NSYN= ",i_sync_step
418 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
421 if (index(controlcard,'TLIST').gt.0) then
423 call card_concat(controlcard1)
424 read(controlcard1,*) (remd_t(i),i=1,nrep)
425 if(me.eq.king.or..not.out1file)
426 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
429 if (index(controlcard,'MLIST').gt.0) then
431 call card_concat(controlcard1)
432 read(controlcard1,*) (remd_m(i),i=1,nrep)
433 if(me.eq.king.or..not.out1file) then
434 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
437 iremd_m_total=iremd_m_total+remd_m(i)
439 write (iout,*) 'Total number of replicas ',iremd_m_total
444 & "Adaptive (PMF-biased) umbrella sampling will be run"
447 if(me.eq.king.or..not.out1file)
448 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
451 c--------------------------------------------------------------------------
452 subroutine read_MDpar
458 include 'COMMON.IOUNITS'
459 include 'COMMON.TIME1'
461 include 'COMMON.QRESTR'
463 include 'COMMON.LANGEVIN'
466 include 'COMMON.LANGEVIN.lang0.5diag'
468 include 'COMMON.LANGEVIN.lang0'
471 include 'COMMON.INTERACT'
472 include 'COMMON.NAMES'
474 include 'COMMON.SETUP'
475 include 'COMMON.CONTROL'
476 include 'COMMON.SPLITELE'
477 include 'COMMON.FFIELD'
479 character*320 controlcard
483 call card_concat(controlcard)
484 call readi(controlcard,"NSTEP",n_timestep,1000000)
485 call readi(controlcard,"NTWE",ntwe,100)
486 call readi(controlcard,"NTWX",ntwx,1000)
487 call readi(controlcard,"REST_FREQ",irest_freq,1000)
488 call reada(controlcard,"DT",d_time,1.0d-1)
489 call reada(controlcard,"DVMAX",dvmax,2.0d1)
490 call reada(controlcard,"DAMAX",damax,1.0d1)
491 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
492 call readi(controlcard,"LANG",lang,0)
493 RESPA = index(controlcard,"RESPA") .gt. 0
494 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
495 ntime_split0=ntime_split
496 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
497 ntime_split0=ntime_split
498 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
499 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
500 rest = index(controlcard,"REST").gt.0
501 tbf = index(controlcard,"TBF").gt.0
502 usampl = index(controlcard,"USAMPL").gt.0
503 scale_umb = index(controlcard,"SCALE_UMB").gt.0
504 adaptive = index(controlcard,"ADAPTIVE").gt.0
505 mdpdb = index(controlcard,"MDPDB").gt.0
506 call reada(controlcard,"T_BATH",t_bath,300.0d0)
507 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
508 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
509 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
510 if (count_reset_moment.eq.0) count_reset_moment=1000000000
511 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
512 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
513 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
514 if (count_reset_vel.eq.0) count_reset_vel=1000000000
515 large = index(controlcard,"LARGE").gt.0
516 print_compon = index(controlcard,"PRINT_COMPON").gt.0
517 rattle = index(controlcard,"RATTLE").gt.0
518 preminim = index(controlcard,"PREMINIM").gt.0
519 if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then
523 write (iout,*) "PREMINIM ",preminim
525 if (index(controlcard,'CART').gt.0) dccart=.true.
526 write (iout,*) "dccart ",dccart
527 write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR')
528 if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim
530 c if performing umbrella sampling, fragments constrained are read from the fragment file
533 write (iout,*) "Umbrella sampling will be run"
534 if (scale_umb.and.adaptive) then
535 write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
536 write (iout,*) "Select one of those and re-run the job."
539 if (scale_umb) write (iout,*)
540 &"Umbrella-restraint force constants will be scaled by temperature"
544 if(me.eq.king.or..not.out1file) then
546 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
548 write (iout,'(a)') "The units are:"
549 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
550 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
551 & " acceleration: angstrom/(48.9 fs)**2"
552 write (iout,'(a)') "energy: kcal/mol, temperature: K"
554 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
555 write (iout,'(a60,f10.5,a)')
556 & "Initial time step of numerical integration:",d_time,
558 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
559 write (iout,'(a60,f10.5,a)') "Cutoff on interactions",r_cut_int,
561 write(iout,'(a60,i5)')"Frequency of updating interaction list",
563 write(iout,'(a60,i5)')"Restart writing frequency",irest_freq
565 write (iout,'(2a,i4,a)')
566 & "A-MTS algorithm used; initial time step for fast-varying",
567 & " short-range forces split into",ntime_split," steps."
568 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
569 & r_cut_respa," lambda",rlamb
571 write (iout,'(2a,f10.5)')
572 & "Maximum acceleration threshold to reduce the time step",
573 & "/increase split number:",damax
574 write (iout,'(2a,f10.5)')
575 & "Maximum predicted energy drift to reduce the timestep",
576 & "/increase split number:",edriftmax
577 write (iout,'(a60,f10.5)')
578 & "Maximum velocity threshold to reduce velocities:",dvmax
579 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
580 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
581 if (rattle) write (iout,'(a60)')
582 & "Rattle algorithm used to constrain the virtual bonds"
586 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
587 call reada(controlcard,"RWAT",rwat,1.4d0)
588 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
589 surfarea=index(controlcard,"SURFAREA").gt.0
590 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
591 if(me.eq.king.or..not.out1file)then
592 write (iout,'(/a,$)') "Langevin dynamics calculation"
595 & " with direct integration of Langevin equations"
596 else if (lang.eq.2) then
597 write (iout,'(a/)') " with TINKER stochasic MD integrator"
598 else if (lang.eq.3) then
599 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
600 else if (lang.eq.4) then
601 write (iout,'(a/)') " in overdamped mode"
603 write (iout,'(//a,i5)')
604 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
607 write (iout,'(a60,f10.5)') "Temperature:",t_bath
608 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
609 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
610 write (iout,'(a60,f10.5)')
611 & "Scaling factor of the friction forces:",scal_fric
612 if (surfarea) write (iout,'(2a,i10,a)')
613 & "Friction coefficients will be scaled by solvent-accessible",
614 & " surface area every",reset_fricmat," steps."
616 c Calculate friction coefficients and bounds of stochastic forces
617 eta=6*pi*cPoise*etawat
618 if(me.eq.king.or..not.out1file)
619 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
621 gamp=scal_fric*(pstok+rwat)*eta
622 stdfp=dsqrt(2*Rb*t_bath/d_time)
624 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
625 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
627 if(me.eq.king.or..not.out1file)then
628 write (iout,'(/2a/)')
629 & "Radii of site types and friction coefficients and std's of",
630 & " stochastic forces of fully exposed sites"
631 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
633 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
634 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
638 if(me.eq.king.or..not.out1file)then
639 write (iout,'(a)') "Berendsen bath calculation"
640 write (iout,'(a60,f10.5)') "Temperature:",t_bath
641 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
643 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
644 & count_reset_moment," steps"
646 & write (iout,'(a,i10,a)')
647 & "Velocities will be reset at random every",count_reset_vel,
651 if(me.eq.king.or..not.out1file)
652 & write (iout,'(a31)') "Microcanonical mode calculation"
654 if(me.eq.king.or..not.out1file)then
655 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
657 write(iout,*) "MD running with constraints."
658 write(iout,*) "Equilibration time ", eq_time, " mtus."
659 write(iout,*) "Constraining ", nfrag," fragments."
660 write(iout,*) "Length of each fragment, weight and q0:"
662 write (iout,*) "Set of restraints #",iset
664 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
665 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
667 write(iout,*) "constraints between ", npair, "fragments."
668 write(iout,*) "constraint pairs, weights and q0:"
670 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
671 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
673 write(iout,*) "angle constraints within ", nfrag_back,
674 & "backbone fragments."
676 write(iout,*) "fragment, weights, q0:"
678 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
679 & ifrag_back(2,i,iset),
680 & wfrag_back(1,i,iset),qin_back(1,i,iset),
681 & wfrag_back(2,i,iset),qin_back(2,i,iset),
682 & wfrag_back(3,i,iset),qin_back(3,i,iset)
685 write(iout,*) "fragment, weights:"
687 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
688 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
689 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
693 iset=mod(kolor,nset)+1
696 if(me.eq.king.or..not.out1file)
697 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
700 c------------------------------------------------------------------------------
703 C Read molecular data.
709 integer error_msg,ierror,ierr,ierrcode
711 include 'COMMON.IOUNITS'
714 include 'COMMON.INTERACT'
715 include 'COMMON.LOCAL'
716 include 'COMMON.NAMES'
717 include 'COMMON.CHAIN'
718 include 'COMMON.FFIELD'
719 include 'COMMON.SBRIDGE'
720 include 'COMMON.HEADER'
721 include 'COMMON.CONTROL'
722 include 'COMMON.SAXS'
723 include 'COMMON.DBASE'
724 include 'COMMON.THREAD'
725 include 'COMMON.CONTACTS'
726 include 'COMMON.TORCNSTR'
727 include 'COMMON.TIME1'
728 include 'COMMON.BOUNDS'
730 include 'COMMON.SETUP'
731 include 'COMMON.SHIELD'
732 character*4 sequence(maxres)
734 double precision x(maxvar)
735 character*256 pdbfile
736 character*400 weightcard
737 character*80 weightcard_t,ucase
738 integer itype_pdb(maxres)
739 common /pizda/ itype_pdb
740 logical seq_comp,fail
741 double precision energia(0:n_ene)
742 double precision secprob(3,maxdih_constr)
744 double precision phihel,phibet,sigmahel,sigmabet
745 integer iti,nsi,maxsi
749 integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp
750 double precision sumv
752 C Read PDB structure if applicable
754 if (indpdb.gt.0 .or. pdbref) then
755 read(inp,'(a)') pdbfile
756 if(me.eq.king.or..not.out1file)
757 & write (iout,'(2a)') 'PDB data will be read from file ',
758 & pdbfile(:ilen(pdbfile))
759 open(ipdbin,file=pdbfile,status='old',err=33)
761 33 write (iout,'(a)') 'Error opening PDB file.'
772 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
773 & (crefjlee(j,i+nres),j=1,3)
776 if(me.eq.king.or..not.out1file)
777 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
778 & ' nstart_sup=',nstart_sup
780 itype_pdb(i)=itype(i)
784 nct=nstart_sup+nsup-1
785 call contact(.false.,ncont_ref,icont_ref,co)
788 if(me.eq.king.or..not.out1file)
789 & write(iout,*)'Adding sidechains'
793 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
796 do while (fail.and.nsi.le.maxsi)
797 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
800 if(fail) write(iout,*)'Adding sidechain failed for res ',
801 & i,' after ',nsi,' trials'
806 if (indpdb.eq.0) then
807 C Read sequence if not taken from the pdb file.
809 c print *,'nres=',nres
810 if (iscode.gt.0) then
811 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
813 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
815 C Convert sequence to numeric code
817 itype(i)=rescode(i,sequence(i),iscode)
821 c print '(20i4)',(itype(i),i=1,nres)
824 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
826 if (itype(i).eq.ntyp1) then
830 else if (iabs(itype(i+1)).ne.20) then
832 else if (iabs(itype(i)).ne.20) then
839 if(me.eq.king.or..not.out1file)then
840 write (iout,*) "ITEL"
842 write (iout,*) i,itype(i),itel(i)
844 c print *,'Call Read_Bridge.'
848 cd print *,'NNT=',NNT,' NCT=',NCT
849 call seq2chains(nres,itype,nchain,chain_length,chain_border,
852 chain_border1(2,1)=chain_border(2,1)+1
854 chain_border1(1,i)=chain_border(1,i)-1
855 chain_border1(2,i)=chain_border(2,i)+1
857 if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
858 chain_border1(2,nchain)=nres
859 write(iout,*) "nres",nres," nchain",nchain
861 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
862 & chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
864 call chain_symmetry(nchain,nres,itype,chain_border,
865 & chain_length,npermchain,tabpermchain,nchain_group,nequiv,
868 c write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
872 write(iout,*) "residue permutations"
874 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
878 if (itype(1).eq.ntyp1) nnt=2
879 if (itype(nres).eq.ntyp1) nct=nct-1
880 write (iout,*) "nnt",nnt," nct",nct
883 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
884 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
886 c print*, 'init_dfa_vars finished!'
888 c print*, 'read_dfa_info finished!'
892 if(me.eq.king.or..not.out1file)
893 & write (iout,'(a,i3)') 'nsup=',nsup
895 if (nsup.le.(nct-nnt+1)) then
896 do i=0,nct-nnt+1-nsup
897 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
903 & 'Error - sequences to be superposed do not match.'
906 do i=0,nsup-(nct-nnt+1)
907 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
909 nstart_sup=nstart_sup+i
915 & 'Error - sequences to be superposed do not match.'
918 if (nsup.eq.0) nsup=nct-nnt
919 if (nstart_sup.eq.0) nstart_sup=nnt
920 if (nstart_seq.eq.0) nstart_seq=nnt
921 if(me.eq.king.or..not.out1file)
922 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
923 & ' nstart_seq=',nstart_seq
926 C 8/13/98 Set limits to generating the dihedral angles
931 read (inp,*) ndih_constr
932 write (iout,*) "ndih_constr",ndih_constr
933 if (ndih_constr.gt.0) then
936 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
939 if(me.eq.king.or..not.out1file)then
942 & 'There are',ndih_constr,' restraints on gamma angles.'
944 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
951 phi0(i)=deg2rad*phi0(i)
952 drange(i)=deg2rad*drange(i)
954 C if(me.eq.king.or..not.out1file)
955 C & write (iout,*) 'FTORS',ftors
958 phibound(1,ii) = phi0(i)-drange(i)
959 phibound(2,ii) = phi0(i)+drange(i)
962 if (me.eq.king .or. .not.out1file) then
964 write (iout,'(a)') 'Boundaries in gamma angle sampling:'
966 write (iout,'(a3,i5,2f10.1)')
967 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
972 else if (ndih_constr.lt.0) then
974 call card_concat(weightcard)
975 call reada(weightcard,"PHIHEL",phihel,50.0D0)
976 call reada(weightcard,"PHIBET",phibet,180.0D0)
977 call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
978 call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
979 call reada(weightcard,"WDIHC",wdihc,0.591D0)
980 write (iout,*) "Weight of dihedral angle restraints",wdihc
981 read(inp,'(9x,3f7.3)')
982 c & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
983 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
984 write (iout,*) "The secprob array"
986 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
991 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
992 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
993 ndih_constr=ndih_constr+1
994 idih_constr(ndih_constr)=i
995 iconstr_dih(i)=ndih_constr
998 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
999 sumv=sumv+vpsipred(j,ndih_constr)
1001 if (sumv.gt.0.0d0) then
1003 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
1006 vpsipred(1,ndih_constr)=1.0d0
1007 vpsipred(2,ndih_constr)=0.0d0
1008 vpsipred(3,ndih_constr)=0.0d0
1010 phibound(1,ndih_constr)=phihel*deg2rad
1011 phibound(2,ndih_constr)=phibet*deg2rad
1012 sdihed(1,ndih_constr)=sigmahel*deg2rad
1013 sdihed(2,ndih_constr)=sigmabet*deg2rad
1017 if(me.eq.king.or..not.out1file)then
1020 & 'There are',ndih_constr,
1021 & ' bimodal restraints on gamma angles.'
1023 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
1024 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
1025 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
1026 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
1027 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
1028 & (vpsipred(j,i),j=1,3)
1035 C first setting the theta boundaries to 0 to pi
1036 C this mean that there is no energy penalty for any angle occuring this can be applied
1037 C for generate random conformation but is not implemented in this way
1040 C thetabound(2,i)=pi
1042 C begin reading theta constrains this is quartic constrains allowing to
1043 C have smooth second derivative
1044 if (with_theta_constr) then
1045 C with_theta_constr is keyword allowing for occurance of theta constrains
1046 read (inp,*) ntheta_constr
1047 C ntheta_constr is the number of theta constrains
1048 if (ntheta_constr.gt.0) then
1049 C read (inp,*) ftors
1050 read (inp,*) (itheta_constr(i),theta_constr0(i),
1051 & theta_drange(i),for_thet_constr(i),
1052 & i=1,ntheta_constr)
1053 C the above code reads from 1 to ntheta_constr
1054 C itheta_constr(i) residue i for which is theta_constr
1055 C theta_constr0 the global minimum value
1056 C theta_drange is range for which there is no energy penalty
1057 C for_thet_constr is the force constant for quartic energy penalty
1059 if(me.eq.king.or..not.out1file)then
1061 & 'There are',ntheta_constr,' constraints on phi angles.'
1062 do i=1,ntheta_constr
1063 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1065 & for_thet_constr(i)
1068 do i=1,ntheta_constr
1069 theta_constr0(i)=deg2rad*theta_constr0(i)
1070 theta_drange(i)=deg2rad*theta_drange(i)
1072 C if(me.eq.king.or..not.out1file)
1073 C & write (iout,*) 'FTORS',ftors
1074 C do i=1,ntheta_constr
1075 C ii = itheta_constr(i)
1076 C thetabound(1,ii) = phi0(i)-drange(i)
1077 C thetabound(2,ii) = phi0(i)+drange(i)
1079 endif ! ntheta_constr.gt.0
1080 endif! with_theta_constr
1082 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1083 C write (iout,*) "with_dihed_constr ",with_dihed_constr
1084 c--- Zscore rms -------
1085 if (nz_start.eq.0) nz_start=nnt
1086 if (nz_end.eq.0 .and. nsup.gt.0) then
1088 else if (nz_end.eq.0) then
1091 if(me.eq.king.or..not.out1file)then
1092 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1093 write (iout,*) 'IZ_SC=',iz_sc
1095 c----------------------
1098 if (.not.pdbref) then
1099 call read_angles(inp,*38)
1102 38 write (iout,'(a)') 'Error reading reference structure.'
1104 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1105 stop 'Error reading reference structure'
1107 39 call chainbuild_extconf
1109 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1118 call contact(.true.,ncont_ref,icont_ref,co)
1122 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1124 if (constr_dist.gt.0) call read_dist_constr
1125 c write (iout,*) "After read_dist_constr nhpb",nhpb
1126 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1128 call NMRpeak_partition
1129 if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1131 if(me.eq.king.or..not.out1file)
1132 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1135 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1137 if(me.eq.king.or..not.out1file)
1138 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1139 & icont_ref(1,i),' ',
1140 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1143 write (iout,*) "calling read_saxs_consrtr",nsaxs
1144 if (nsaxs.gt.0) call read_saxs_constr
1145 c write (iout,*) "After read_saxs_constr"
1147 if (constr_homology.gt.0) then
1148 c write (iout,*) "Calling read_constr_homology"
1150 call read_constr_homology
1151 if (indpdb.gt.0 .or. pdbref) then
1154 c(j,i)=crefjlee(j,i)
1155 cref(j,i)=crefjlee(j,i)
1160 write (iout,*) "sc_loc_geom: Array C"
1162 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1163 & (c(j,i+nres),j=1,3)
1165 write (iout,*) "Array Cref"
1167 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1168 & (cref(j,i+nres),j=1,3)
1171 call int_from_cart1(.false.)
1172 call sc_loc_geom(.false.)
1174 thetaref(i)=theta(i)
1179 dc(j,i)=c(j,i+1)-c(j,i)
1180 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1185 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1186 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1191 if (start_from_model) then
1194 read(inp,'(a)',end=332,err=332) pdbfile
1195 if (me.eq.king .or. .not. out1file)
1196 & write (iout,'(a,5x,a)') 'Opening PDB file',
1197 & pdbfile(:ilen(pdbfile))
1198 open(ipdbin,file=pdbfile,status='old',err=336)
1200 336 write (iout,'(a,5x,a)') 'Error opening PDB file',
1201 & pdbfile(:ilen(pdbfile))
1208 call readpdb_template(nmodel_start+1)
1210 if (nres.ge.nres_temp) then
1211 nmodel_start=nmodel_start+1
1212 pdbfiles_chomo(nmodel_start)=pdbfile
1215 chomo(j,i,nmodel_start)=c(j,i)
1219 if (me.eq.king .or. .not. out1file)
1220 & write (iout,'(a,2i7,1x,a)')
1221 & "Different number of residues",nres_temp,nres,
1227 if (nmodel_start.eq.0) then
1228 if (me.eq.king .or. .not. out1file)
1229 & write (iout,'(a)')
1230 & "No valid starting model found START_FROM_MODELS is OFF"
1231 start_from_model=.false.
1233 write (iout,*) "nmodel_start",nmodel_start
1239 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1240 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1241 & modecalc.ne.10) then
1242 C If input structure hasn't been supplied from the PDB file read or generate
1244 if (iranconf.eq.0 .and. .not. extconf .and. .not.
1245 & start_from_model) then
1246 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1247 & write (iout,'(a)') 'Initial geometry will be read in.'
1249 read(inp,'(8f10.5)',end=36,err=36)
1250 & ((c(l,k),l=1,3),k=1,nres),
1251 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1252 if (nnt.gt.1) c(:,nres+1)=c(:,1)
1253 if (nct.lt.nres) c(:,2*nres)=c(:,nres)
1254 c write (iout,*) "Exit READ_CART"
1255 c write (iout,'(8f10.5)')
1256 c & ((c(l,k),l=1,3),k=1,nres),
1257 c & ((c(l,k+nres),l=1,3),k=nnt,nct)
1264 dc(j,i)=c(j,i+1)-c(j,i)
1265 c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1269 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1271 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1272 c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1277 c dc_norm(j,i+nres)=0.0d0
1281 call int_from_cart1(.true.)
1282 write (iout,*) "Finish INT_TO_CART"
1283 c write (iout,*) "DC"
1285 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1286 c & (dc(j,i+nres),j=1,3)
1292 call read_angles(inp,*36)
1294 call chainbuild_extconf
1297 36 write (iout,'(a)') 'Error reading angle file.'
1299 call mpi_finalize( MPI_COMM_WORLD,IERR )
1301 stop 'Error reading angle file.'
1303 else if (extconf) then
1304 if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1305 & write (iout,'(a)') 'Extended chain initial geometry.'
1307 theta(i)=90d0*deg2rad
1310 phi(i)=180d0*deg2rad
1313 alph(i)=110d0*deg2rad
1316 omeg(i)=-120d0*deg2rad
1317 if (itype(i).le.0) omeg(i)=-omeg(i)
1320 call chainbuild_extconf
1321 else if (.not. start_from_model) then
1322 if(me.eq.king.or..not.out1file)
1323 & write (iout,'(a)') 'Random-generated initial geometry.'
1326 if (me.eq.king .or. fg_rank.eq.0 .and. (
1327 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1331 call gen_rand_conf(itmp,*30)
1333 30 write (iout,*) 'Failed to generate random conformation',
1334 & ', itrial=',itrial
1335 write (*,*) 'Processor:',me,
1336 & ' Failed to generate random conformation',
1346 write (iout,'(a,i3,a)') 'Processor:',me,
1347 & ' error in generating random conformation.'
1348 write (*,'(a,i3,a)') 'Processor:',me,
1349 & ' error in generating random conformation.'
1352 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1357 & ' error in generating random conformation.'
1362 elseif (modecalc.eq.4) then
1363 read (inp,'(a)') intinname
1364 open (intin,file=intinname,status='old',err=333)
1365 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1366 & write (iout,'(a)') 'intinname',intinname
1367 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1369 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1371 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1373 stop 'Error opening angle file.'
1377 C Generate distance constraints, if the PDB structure is to be regularized.
1378 if (nthread.gt.0) then
1379 call read_threadbase
1382 if (me.eq.king .or. .not. out1file)
1384 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1385 write (iout,'(/a,i3,a)')
1386 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1387 write (iout,'(20i4)') (iss(i),i=1,ns)
1389 write(iout,*)"Running with dynamic disulfide-bond formation"
1391 write (iout,'(/a/)') 'Pre-formed links are:'
1397 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1398 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1404 if (ns.gt.0.and.dyn_ss) then
1408 forcon(i-nss)=forcon(i)
1415 dyn_ss_mask(iss(i))=.true.
1420 c write (iout,*) "DC"
1422 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1423 c & (dc(j,i+nres),j=1,3)
1425 if (i2ndstr.gt.0) call secstrp2dihc
1426 c call geom_to_var(nvar,x)
1427 c call etotal(energia(0))
1428 c call enerprint(energia(0))
1429 c call briefout(0,etot)
1431 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1432 cd write (iout,'(a)') 'Variable list:'
1433 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1435 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1436 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1437 & 'Processor',myrank,': end reading molecular data.'
1442 c--------------------------------------------------------------------------
1443 logical function seq_comp(itypea,itypeb,length)
1445 integer length,itypea(length),itypeb(length)
1448 if (itypea(i).ne.itypeb(i)) then
1456 c-----------------------------------------------------------------------------
1457 subroutine read_bridge
1458 C Read information about disulfide bridges.
1460 include 'DIMENSIONS'
1465 include 'COMMON.IOUNITS'
1466 include 'COMMON.GEO'
1467 include 'COMMON.VAR'
1468 include 'COMMON.INTERACT'
1469 include 'COMMON.LOCAL'
1470 include 'COMMON.NAMES'
1471 include 'COMMON.CHAIN'
1472 include 'COMMON.FFIELD'
1473 include 'COMMON.SBRIDGE'
1474 include 'COMMON.HEADER'
1475 include 'COMMON.CONTROL'
1476 include 'COMMON.DBASE'
1477 include 'COMMON.THREAD'
1478 include 'COMMON.TIME1'
1479 include 'COMMON.SETUP'
1481 C Read bridging residues.
1482 read (inp,*) ns,(iss(i),i=1,ns)
1483 c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
1490 if(me.eq.king.or..not.out1file)
1491 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1492 C Check whether the specified bridging residues are cystines.
1494 if (itype(iss(i)).ne.1) then
1495 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1496 & 'Do you REALLY think that the residue ',
1497 & restyp(itype(iss(i))),i,
1498 & ' can form a disulfide bridge?!!!'
1499 write (*,'(2a,i3,a)')
1500 & 'Do you REALLY think that the residue ',
1501 & restyp(itype(iss(i))),i,
1502 & ' can form a disulfide bridge?!!!'
1504 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1509 C Read preformed bridges.
1511 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1513 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1516 C Check if the residues involved in bridges are in the specified list of
1517 C bridging residues.
1520 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1521 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1522 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1523 & ' contains residues present in other pairs.'
1524 write (*,'(a,i3,a)') 'Disulfide pair',i,
1525 & ' contains residues present in other pairs.'
1527 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1533 if (ihpb(i).eq.iss(j)) goto 10
1535 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1538 if (jhpb(i).eq.iss(j)) goto 20
1540 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1546 ihpb(i)=ihpb(i)+nres
1547 jhpb(i)=jhpb(i)+nres
1553 c----------------------------------------------------------------------------
1554 subroutine read_x(kanal,*)
1556 include 'DIMENSIONS'
1557 include 'COMMON.GEO'
1558 include 'COMMON.VAR'
1559 include 'COMMON.CHAIN'
1560 include 'COMMON.IOUNITS'
1561 include 'COMMON.CONTROL'
1562 include 'COMMON.LOCAL'
1563 include 'COMMON.INTERACT'
1564 integer i,j,k,l,kanal
1565 c Read coordinates from input
1567 read(kanal,'(8f10.5)',end=10,err=10)
1568 & ((c(l,k),l=1,3),k=1,nres),
1569 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1572 c(j,2*nres)=c(j,nres)
1574 call int_from_cart1(.false.)
1577 dc(j,i)=c(j,i+1)-c(j,i)
1578 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1582 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1584 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1585 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1593 c----------------------------------------------------------------------------
1594 subroutine read_threadbase
1596 include 'DIMENSIONS'
1597 include 'COMMON.IOUNITS'
1598 include 'COMMON.GEO'
1599 include 'COMMON.VAR'
1600 include 'COMMON.INTERACT'
1601 include 'COMMON.LOCAL'
1602 include 'COMMON.NAMES'
1603 include 'COMMON.CHAIN'
1604 include 'COMMON.FFIELD'
1605 include 'COMMON.SBRIDGE'
1606 include 'COMMON.HEADER'
1607 include 'COMMON.CONTROL'
1608 include 'COMMON.DBASE'
1609 include 'COMMON.THREAD'
1610 include 'COMMON.TIME1'
1612 double precision dist
1613 C Read pattern database for threading.
1614 read (icbase,*) nseq
1616 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1617 & nres_base(2,i),nres_base(3,i)
1618 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1620 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1621 c & nres_base(2,i),nres_base(3,i)
1622 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1626 if (weidis.eq.0.0D0) weidis=0.1D0
1635 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1636 write (iout,'(a,i5)') 'nexcl: ',nexcl
1637 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1640 c------------------------------------------------------------------------------
1641 subroutine setup_var
1643 include 'DIMENSIONS'
1644 include 'COMMON.IOUNITS'
1645 include 'COMMON.GEO'
1646 include 'COMMON.VAR'
1647 include 'COMMON.INTERACT'
1648 include 'COMMON.LOCAL'
1649 include 'COMMON.NAMES'
1650 include 'COMMON.CHAIN'
1651 include 'COMMON.FFIELD'
1652 include 'COMMON.SBRIDGE'
1653 include 'COMMON.HEADER'
1654 include 'COMMON.CONTROL'
1655 include 'COMMON.DBASE'
1656 include 'COMMON.THREAD'
1657 include 'COMMON.TIME1'
1659 C Set up variable list.
1665 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1667 ialph(i,1)=nvar+nside
1671 if (indphi.gt.0) then
1673 else if (indback.gt.0) then
1680 c----------------------------------------------------------------------------
1681 subroutine gen_dist_constr
1682 C Generate CA distance constraints.
1684 include 'DIMENSIONS'
1685 include 'COMMON.IOUNITS'
1686 include 'COMMON.GEO'
1687 include 'COMMON.VAR'
1688 include 'COMMON.INTERACT'
1689 include 'COMMON.LOCAL'
1690 include 'COMMON.NAMES'
1691 include 'COMMON.CHAIN'
1692 include 'COMMON.FFIELD'
1693 include 'COMMON.SBRIDGE'
1694 include 'COMMON.HEADER'
1695 include 'COMMON.CONTROL'
1696 include 'COMMON.DBASE'
1697 include 'COMMON.THREAD'
1698 include 'COMMON.SPLITELE'
1699 include 'COMMON.TIME1'
1700 integer i,j,itype_pdb(maxres)
1701 common /pizda/ itype_pdb
1703 double precision dist
1705 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1706 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1707 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1709 do i=nstart_sup,nstart_sup+nsup-1
1710 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1711 cd & ' seq_pdb', restyp(itype_pdb(i))
1712 do j=i+2,nstart_sup+nsup-1
1713 c 5/24/2020 Adam: Cutoff included to reduce array size
1715 if (dd.gt.r_cut_int) cycle
1717 ihpb(nhpb)=i+nstart_seq-nstart_sup
1718 jhpb(nhpb)=j+nstart_seq-nstart_sup
1723 cd write (iout,'(a)') 'Distance constraints:'
1728 cd if (ii.gt.nres) then
1733 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1734 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1735 cd & dhpb(i),forcon(i)
1739 c----------------------------------------------------------------------------
1742 include 'DIMENSIONS'
1743 include 'COMMON.MAP'
1744 include 'COMMON.IOUNITS'
1746 character*3 angid(4) /'THE','PHI','ALP','OME'/
1747 character*80 mapcard,ucase
1749 read (inp,'(a)') mapcard
1750 mapcard=ucase(mapcard)
1751 if (index(mapcard,'PHI').gt.0) then
1753 else if (index(mapcard,'THE').gt.0) then
1755 else if (index(mapcard,'ALP').gt.0) then
1757 else if (index(mapcard,'OME').gt.0) then
1760 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1761 stop 'Error - illegal variable spec in MAP card.'
1763 call readi (mapcard,'RES1',res1(imap),0)
1764 call readi (mapcard,'RES2',res2(imap),0)
1765 if (res1(imap).eq.0) then
1766 res1(imap)=res2(imap)
1767 else if (res2(imap).eq.0) then
1768 res2(imap)=res1(imap)
1770 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1772 & 'Error - illegal definition of variable group in MAP.'
1773 stop 'Error - illegal definition of variable group in MAP.'
1775 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1776 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1777 call readi(mapcard,'NSTEP',nstep(imap),0)
1778 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1780 & 'Illegal boundary and/or step size specification in MAP.'
1781 stop 'Illegal boundary and/or step size specification in MAP.'
1786 c----------------------------------------------------------------------------
1789 include 'DIMENSIONS'
1790 include 'COMMON.IOUNITS'
1791 include 'COMMON.GEO'
1792 include 'COMMON.CSA'
1793 include 'COMMON.BANK'
1794 include 'COMMON.CONTROL'
1796 character*620 mcmcard
1797 call card_concat(mcmcard)
1799 call readi(mcmcard,'NCONF',nconf,50)
1800 call readi(mcmcard,'NADD',nadd,0)
1801 call readi(mcmcard,'JSTART',jstart,1)
1802 call readi(mcmcard,'JEND',jend,1)
1803 call readi(mcmcard,'NSTMAX',nstmax,500000)
1804 call readi(mcmcard,'N0',n0,1)
1805 call readi(mcmcard,'N1',n1,6)
1806 call readi(mcmcard,'N2',n2,4)
1807 call readi(mcmcard,'N3',n3,0)
1808 call readi(mcmcard,'N4',n4,0)
1809 call readi(mcmcard,'N5',n5,0)
1810 call readi(mcmcard,'N6',n6,10)
1811 call readi(mcmcard,'N7',n7,0)
1812 call readi(mcmcard,'N8',n8,0)
1813 call readi(mcmcard,'N9',n9,0)
1814 call readi(mcmcard,'N14',n14,0)
1815 call readi(mcmcard,'N15',n15,0)
1816 call readi(mcmcard,'N16',n16,0)
1817 call readi(mcmcard,'N17',n17,0)
1818 call readi(mcmcard,'N18',n18,0)
1820 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1822 call readi(mcmcard,'NDIFF',ndiff,2)
1823 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1824 call readi(mcmcard,'IS1',is1,1)
1825 call readi(mcmcard,'IS2',is2,8)
1826 call readi(mcmcard,'NRAN0',nran0,4)
1827 call readi(mcmcard,'NRAN1',nran1,2)
1828 call readi(mcmcard,'IRR',irr,1)
1829 call readi(mcmcard,'NSEED',nseed,20)
1830 call readi(mcmcard,'NTOTAL',ntotal,10000)
1831 call reada(mcmcard,'CUT1',cut1,2.0d0)
1832 call reada(mcmcard,'CUT2',cut2,5.0d0)
1833 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1834 call readi(mcmcard,'ICMAX',icmax,3)
1835 call readi(mcmcard,'IRESTART',irestart,0)
1836 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1839 call reada(mcmcard,'DELE',dele,20.0d0)
1840 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1841 call readi(mcmcard,'IREF',iref,0)
1842 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1843 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1844 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1845 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1846 write (iout,*) "NCONF_IN",nconf_in
1849 c----------------------------------------------------------------------------
1852 include 'DIMENSIONS'
1853 include 'COMMON.MCM'
1854 include 'COMMON.MCE'
1855 include 'COMMON.IOUNITS'
1857 character*320 mcmcard
1859 call card_concat(mcmcard)
1860 call readi(mcmcard,'MAXACC',maxacc,100)
1861 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1862 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1863 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1864 call readi(mcmcard,'MAXREPM',maxrepm,200)
1865 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1866 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1867 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1868 call reada(mcmcard,'E_UP',e_up,5.0D0)
1869 call reada(mcmcard,'DELTE',delte,0.1D0)
1870 call readi(mcmcard,'NSWEEP',nsweep,5)
1871 call readi(mcmcard,'NSTEPH',nsteph,0)
1872 call readi(mcmcard,'NSTEPC',nstepc,0)
1873 call reada(mcmcard,'TMIN',tmin,298.0D0)
1874 call reada(mcmcard,'TMAX',tmax,298.0D0)
1875 call readi(mcmcard,'NWINDOW',nwindow,0)
1876 call readi(mcmcard,'PRINT_MC',print_mc,0)
1877 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1878 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1879 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1880 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1881 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1882 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1883 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1884 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1885 if (nwindow.gt.0) then
1886 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1888 winlen(i)=winend(i)-winstart(i)+1
1891 if (tmax.lt.tmin) tmax=tmin
1892 if (tmax.eq.tmin) then
1896 if (nstepc.gt.0 .and. nsteph.gt.0) then
1897 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1898 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1900 C Probabilities of different move types
1901 sumpro_type(0)=0.0D0
1902 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1903 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1904 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1905 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1906 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1907 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1908 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1910 print *,'i',i,' sumprotype',sumpro_type(i)
1911 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1912 print *,'i',i,' sumprotype',sumpro_type(i)
1916 c----------------------------------------------------------------------------
1917 subroutine read_minim
1919 include 'DIMENSIONS'
1920 include 'COMMON.MINIM'
1921 include 'COMMON.IOUNITS'
1923 character*320 minimcard
1924 call card_concat(minimcard)
1925 call readi(minimcard,'MAXMIN',maxmin,2000)
1926 call readi(minimcard,'MAXFUN',maxfun,5000)
1927 call readi(minimcard,'MINMIN',minmin,maxmin)
1928 call readi(minimcard,'MINFUN',minfun,maxmin)
1929 call reada(minimcard,'TOLF',tolf,1.0D-2)
1930 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1931 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1932 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1933 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1934 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1935 & 'Options in energy minimization:'
1936 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1937 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1938 & 'MinMin:',MinMin,' MinFun:',MinFun,
1939 & ' TolF:',TolF,' RTolF:',RTolF
1942 c----------------------------------------------------------------------------
1943 subroutine read_angles(kanal,*)
1945 include 'DIMENSIONS'
1946 include 'COMMON.GEO'
1947 include 'COMMON.VAR'
1948 include 'COMMON.CHAIN'
1949 include 'COMMON.IOUNITS'
1950 include 'COMMON.CONTROL'
1952 c Read angles from input
1954 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1955 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1956 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1957 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1960 c 9/7/01 avoid 180 deg valence angle
1961 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1963 theta(i)=deg2rad*theta(i)
1964 phi(i)=deg2rad*phi(i)
1965 alph(i)=deg2rad*alph(i)
1966 omeg(i)=deg2rad*omeg(i)
1971 c----------------------------------------------------------------------------
1972 subroutine reada(rekord,lancuch,wartosc,default)
1974 character*(*) rekord,lancuch
1975 double precision wartosc,default
1978 iread=index(rekord,lancuch)
1979 if (iread.eq.0) then
1983 iread=iread+ilen(lancuch)+1
1984 read (rekord(iread:),*,err=10,end=10) wartosc
1989 c----------------------------------------------------------------------------
1990 subroutine readi(rekord,lancuch,wartosc,default)
1992 character*(*) rekord,lancuch
1993 integer wartosc,default
1996 iread=index(rekord,lancuch)
1997 if (iread.eq.0) then
2001 iread=iread+ilen(lancuch)+1
2002 read (rekord(iread:),*,err=10,end=10) wartosc
2007 c----------------------------------------------------------------------------
2008 subroutine multreadi(rekord,lancuch,tablica,dim,default)
2011 integer tablica(dim),default
2012 character*(*) rekord,lancuch
2019 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2020 if (iread.eq.0) return
2021 iread=iread+ilen(lancuch)+1
2022 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2025 c----------------------------------------------------------------------------
2026 subroutine multreada(rekord,lancuch,tablica,dim,default)
2029 double precision tablica(dim),default
2030 character*(*) rekord,lancuch
2037 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2038 if (iread.eq.0) return
2039 iread=iread+ilen(lancuch)+1
2040 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2043 c----------------------------------------------------------------------------
2044 subroutine openunits
2046 include 'DIMENSIONS'
2050 character*16 form,nodename
2053 include 'COMMON.SETUP'
2054 include 'COMMON.IOUNITS'
2056 include 'COMMON.CONTROL'
2057 integer lenpre,lenpot,ilen,lentmp,npos
2059 character*3 out1file_text,ucase
2064 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2065 call getenv_loc("PREFIX",prefix)
2067 call getenv_loc("POT",pot)
2068 call getenv_loc("DIRTMP",tmpdir)
2069 call getenv_loc("CURDIR",curdir)
2070 call getenv_loc("OUT1FILE",out1file_text)
2071 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2072 out1file_text=ucase(out1file_text)
2073 if (out1file_text(1:1).eq."Y") then
2076 out1file=fg_rank.gt.0
2081 if (lentmp.gt.0) then
2082 write (*,'(80(1h!))')
2083 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2084 write (*,'(80(1h!))')
2085 write (*,*)"All output files will be on node /tmp directory."
2087 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2088 if (me.eq.king) then
2089 write (*,*) "The master node is ",nodename
2090 else if (fg_rank.eq.0) then
2091 write (*,*) "I am the CG slave node ",nodename
2093 write (*,*) "I am the FG slave node ",nodename
2096 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2097 lenpre = lentmp+lenpre+1
2099 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2100 C Get the names and open the input files
2101 #if defined(WINIFL) || defined(WINPGI)
2102 open(1,file=pref_orig(:ilen(pref_orig))//
2103 & '.inp',status='old',readonly,shared)
2104 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2105 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2106 C Get parameter filenames and open the parameter files.
2107 call getenv_loc('BONDPAR',bondname)
2108 open (ibond,file=bondname,status='old',readonly,shared)
2109 call getenv_loc('THETPAR',thetname)
2110 open (ithep,file=thetname,status='old',readonly,shared)
2111 call getenv_loc('ROTPAR',rotname)
2112 open (irotam,file=rotname,status='old',readonly,shared)
2113 call getenv_loc('TORPAR',torname)
2114 open (itorp,file=torname,status='old',readonly,shared)
2115 call getenv_loc('TORDPAR',tordname)
2116 open (itordp,file=tordname,status='old',readonly,shared)
2117 call getenv_loc('FOURIER',fouriername)
2118 open (ifourier,file=fouriername,status='old',readonly,shared)
2119 call getenv_loc('ELEPAR',elename)
2120 open (ielep,file=elename,status='old',readonly,shared)
2121 call getenv_loc('SIDEPAR',sidename)
2122 open (isidep,file=sidename,status='old',readonly,shared)
2123 call getenv_loc('LIPTRANPAR',liptranname)
2124 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2125 #elif (defined CRAY) || (defined AIX)
2126 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2128 c print *,"Processor",myrank," opened file 1"
2129 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2130 c print *,"Processor",myrank," opened file 9"
2131 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2132 C Get parameter filenames and open the parameter files.
2133 call getenv_loc('BONDPAR',bondname)
2134 open (ibond,file=bondname,status='old',action='read')
2135 c print *,"Processor",myrank," opened file IBOND"
2136 call getenv_loc('THETPAR',thetname)
2137 open (ithep,file=thetname,status='old',action='read')
2139 call getenv_loc('THETPARPDB',thetname_pdb)
2140 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2142 c print *,"Processor",myrank," opened file ITHEP"
2143 call getenv_loc('ROTPAR',rotname)
2144 open (irotam,file=rotname,status='old',action='read')
2146 call getenv_loc('ROTPARPDB',rotname_pdb)
2147 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2149 c print *,"Processor",myrank," opened file IROTAM"
2150 call getenv_loc('TORPAR',torname)
2151 open (itorp,file=torname,status='old',action='read')
2152 c print *,"Processor",myrank," opened file ITORP"
2153 call getenv_loc('TORDPAR',tordname)
2154 open (itordp,file=tordname,status='old',action='read')
2155 c print *,"Processor",myrank," opened file ITORDP"
2156 call getenv_loc('SCCORPAR',sccorname)
2157 open (isccor,file=sccorname,status='old',action='read')
2158 c print *,"Processor",myrank," opened file ISCCOR"
2159 call getenv_loc('FOURIER',fouriername)
2160 open (ifourier,file=fouriername,status='old',action='read')
2161 c print *,"Processor",myrank," opened file IFOURIER"
2162 call getenv_loc('ELEPAR',elename)
2163 open (ielep,file=elename,status='old',action='read')
2164 c print *,"Processor",myrank," opened file IELEP"
2165 call getenv_loc('SIDEPAR',sidename)
2166 open (isidep,file=sidename,status='old',action='read')
2167 call getenv_loc('LIPTRANPAR',liptranname)
2168 open (iliptranpar,file=liptranname,status='old',action='read')
2169 c print *,"Processor",myrank," opened file ISIDEP"
2170 c print *,"Processor",myrank," opened parameter files"
2172 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2173 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2174 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2175 C Get parameter filenames and open the parameter files.
2176 call getenv_loc('BONDPAR',bondname)
2177 open (ibond,file=bondname,status='old')
2178 call getenv_loc('THETPAR',thetname)
2179 open (ithep,file=thetname,status='old')
2181 call getenv_loc('THETPARPDB',thetname_pdb)
2182 open (ithep_pdb,file=thetname_pdb,status='old')
2184 call getenv_loc('ROTPAR',rotname)
2185 open (irotam,file=rotname,status='old')
2187 call getenv_loc('ROTPARPDB',rotname_pdb)
2188 open (irotam_pdb,file=rotname_pdb,status='old')
2190 call getenv_loc('TORPAR',torname)
2191 open (itorp,file=torname,status='old')
2192 call getenv_loc('TORDPAR',tordname)
2193 open (itordp,file=tordname,status='old')
2194 call getenv_loc('SCCORPAR',sccorname)
2195 open (isccor,file=sccorname,status='old')
2196 call getenv_loc('FOURIER',fouriername)
2197 open (ifourier,file=fouriername,status='old')
2198 call getenv_loc('ELEPAR',elename)
2199 open (ielep,file=elename,status='old')
2200 call getenv_loc('SIDEPAR',sidename)
2201 open (isidep,file=sidename,status='old')
2202 call getenv_loc('LIPTRANPAR',liptranname)
2203 open (iliptranpar,file=liptranname,status='old')
2205 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2207 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2208 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2209 C Get parameter filenames and open the parameter files.
2210 call getenv_loc('BONDPAR',bondname)
2211 open (ibond,file=bondname,status='old',readonly)
2212 call getenv_loc('THETPAR',thetname)
2213 open (ithep,file=thetname,status='old',readonly)
2214 call getenv_loc('ROTPAR',rotname)
2215 open (irotam,file=rotname,status='old',readonly)
2216 call getenv_loc('TORPAR',torname)
2217 open (itorp,file=torname,status='old',readonly)
2218 call getenv_loc('TORDPAR',tordname)
2219 open (itordp,file=tordname,status='old',readonly)
2220 call getenv_loc('SCCORPAR',sccorname)
2221 open (isccor,file=sccorname,status='old',readonly)
2223 call getenv_loc('THETPARPDB',thetname_pdb)
2224 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2226 call getenv_loc('FOURIER',fouriername)
2227 open (ifourier,file=fouriername,status='old',readonly)
2228 call getenv_loc('ELEPAR',elename)
2229 open (ielep,file=elename,status='old',readonly)
2230 call getenv_loc('SIDEPAR',sidename)
2231 open (isidep,file=sidename,status='old',readonly)
2232 call getenv_loc('LIPTRANPAR',liptranname)
2233 open (iliptranpar,file=liptranname,status='old',action='read')
2235 call getenv_loc('ROTPARPDB',rotname_pdb)
2236 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2240 call getenv_loc('TUBEPAR',tubename)
2241 #if defined(WINIFL) || defined(WINPGI)
2242 open (itube,file=tubename,status='old',readonly,shared)
2243 #elif (defined CRAY) || (defined AIX)
2244 open (itube,file=tubename,status='old',action='read')
2246 open (itube,file=tubename,status='old')
2248 open (itube,file=tubename,status='old',readonly)
2253 C 8/9/01 In the newest version SCp interaction constants are read from a file
2254 C Use -DOLDSCP to use hard-coded constants instead.
2256 call getenv_loc('SCPPAR',scpname)
2257 #if defined(WINIFL) || defined(WINPGI)
2258 open (iscpp,file=scpname,status='old',readonly,shared)
2259 #elif (defined CRAY) || (defined AIX)
2260 open (iscpp,file=scpname,status='old',action='read')
2262 open (iscpp,file=scpname,status='old')
2264 open (iscpp,file=scpname,status='old',readonly)
2267 call getenv_loc('PATTERN',patname)
2268 #if defined(WINIFL) || defined(WINPGI)
2269 open (icbase,file=patname,status='old',readonly,shared)
2270 #elif (defined CRAY) || (defined AIX)
2271 open (icbase,file=patname,status='old',action='read')
2273 open (icbase,file=patname,status='old')
2275 open (icbase,file=patname,status='old',readonly)
2278 C Open output file only for CG processes
2279 c print *,"Processor",myrank," fg_rank",fg_rank
2280 if (fg_rank.eq.0) then
2282 if (nodes.eq.1) then
2285 npos = dlog10(dfloat(nodes-1))+1
2287 if (npos.lt.3) npos=3
2288 write (liczba,'(i1)') npos
2289 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2291 write (liczba,form) me
2292 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2293 & liczba(:ilen(liczba))
2294 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2296 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2298 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2299 & liczba(:ilen(liczba))//'.mol2'
2300 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2301 & liczba(:ilen(liczba))//'.stat'
2303 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2304 & //liczba(:ilen(liczba))//'.stat')
2305 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2308 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2309 & liczba(:ilen(liczba))//'.const'
2314 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2315 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2316 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2317 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2318 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2320 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2322 rest2name=prefix(:ilen(prefix))//'.rst'
2324 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2327 #if defined(AIX) || defined(PGI) || defined(CRAY)
2328 if (me.eq.king .or. .not. out1file) then
2329 open(iout,file=outname,status='unknown')
2331 open(iout,file="/dev/null",status="unknown")
2335 if (fg_rank.gt.0) then
2336 write (liczba,'(i3.3)') myrank/nfgtasks
2337 write (ll,'(bz,i3.3)') fg_rank
2338 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2344 open(igeom,file=intname,status='unknown',position='append')
2345 open(ipdb,file=pdbname,status='unknown')
2346 open(imol2,file=mol2name,status='unknown')
2347 open(istat,file=statname,status='unknown',position='append')
2349 c1out open(iout,file=outname,status='unknown')
2352 if (me.eq.king .or. .not.out1file)
2353 & open(iout,file=outname,status='unknown')
2356 if (fg_rank.gt.0) then
2357 write (liczba,'(i3.3)') myrank/nfgtasks
2358 write (ll,'(bz,i3.3)') fg_rank
2359 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2365 open(igeom,file=intname,status='unknown',access='append')
2366 open(ipdb,file=pdbname,status='unknown')
2367 open(imol2,file=mol2name,status='unknown')
2368 open(istat,file=statname,status='unknown',access='append')
2370 c1out open(iout,file=outname,status='unknown')
2373 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2374 csa_seed=prefix(:lenpre)//'.CSA.seed'
2375 csa_history=prefix(:lenpre)//'.CSA.history'
2376 csa_bank=prefix(:lenpre)//'.CSA.bank'
2377 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2378 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2379 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2380 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2381 csa_int=prefix(:lenpre)//'.int'
2382 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2383 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2384 csa_in=prefix(:lenpre)//'.CSA.in'
2385 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2388 write (iout,'(80(1h-))')
2389 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2390 write (iout,'(80(1h-))')
2391 write (iout,*) "Input file : ",
2392 & pref_orig(:ilen(pref_orig))//'.inp'
2393 write (iout,*) "Output file : ",
2394 & outname(:ilen(outname))
2396 write (iout,*) "Sidechain potential file : ",
2397 & sidename(:ilen(sidename))
2399 write (iout,*) "SCp potential file : ",
2400 & scpname(:ilen(scpname))
2402 write (iout,*) "Electrostatic potential file : ",
2403 & elename(:ilen(elename))
2404 write (iout,*) "Cumulant coefficient file : ",
2405 & fouriername(:ilen(fouriername))
2406 write (iout,*) "Torsional parameter file : ",
2407 & torname(:ilen(torname))
2408 write (iout,*) "Double torsional parameter file : ",
2409 & tordname(:ilen(tordname))
2410 write (iout,*) "SCCOR parameter file : ",
2411 & sccorname(:ilen(sccorname))
2412 write (iout,*) "Bond & inertia constant file : ",
2413 & bondname(:ilen(bondname))
2414 write (iout,*) "Bending-torsion parameter file : ",
2415 & thetname(:ilen(thetname))
2416 write (iout,*) "Rotamer parameter file : ",
2417 & rotname(:ilen(rotname))
2418 write (iout,*) "Thetpdb parameter file : ",
2419 & thetname_pdb(:ilen(thetname_pdb))
2420 write (iout,*) "Threading database : ",
2421 & patname(:ilen(patname))
2423 &write (iout,*)" DIRTMP : ",
2425 write (iout,'(80(1h-))')
2429 c----------------------------------------------------------------------------
2430 subroutine card_concat(card)
2431 implicit real*8 (a-h,o-z)
2432 include 'DIMENSIONS'
2433 include 'COMMON.IOUNITS'
2435 character*80 karta,ucase
2437 read (inp,'(a)') karta
2440 do while (karta(80:80).eq.'&')
2441 card=card(:ilen(card)+1)//karta(:79)
2442 read (inp,'(a)') karta
2445 card=card(:ilen(card)+1)//karta
2448 c------------------------------------------------------------------------------
2451 include 'DIMENSIONS'
2452 include 'COMMON.CHAIN'
2453 include 'COMMON.IOUNITS'
2454 include 'COMMON.CONTROL'
2456 include 'COMMON.QRESTR'
2458 include 'COMMON.LAGRANGE.5diag'
2460 include 'COMMON.LAGRANGE'
2463 open(irest2,file=rest2name,status='unknown')
2464 read(irest2,*) totT,EK,potE,totE,t_bath
2467 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2470 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2473 read (irest2,*) iset
2478 c------------------------------------------------------------------------------
2479 subroutine read_fragments
2481 include 'DIMENSIONS'
2485 include 'COMMON.SETUP'
2486 include 'COMMON.CHAIN'
2487 include 'COMMON.IOUNITS'
2489 include 'COMMON.QRESTR'
2490 include 'COMMON.CONTROL'
2492 read(inp,*) nset,nfrag,npair,nfrag_back
2493 loc_qlike=(nfrag_back.lt.0)
2494 nfrag_back=iabs(nfrag_back)
2495 if(me.eq.king.or..not.out1file)
2496 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2497 & " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2499 read(inp,*) mset(iset)
2501 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2503 if(me.eq.king.or..not.out1file)
2504 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2505 & ifrag(2,i,iset), qinfrag(i,iset)
2508 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2510 if(me.eq.king.or..not.out1file)
2511 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2512 & ipair(2,i,iset), qinpair(i,iset)
2516 read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2517 & wfrag_back(2,i,iset),qin_back(2,i,iset),
2518 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2519 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2520 if(me.eq.king.or..not.out1file)
2521 & write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2522 & wfrag_back(2,i,iset),qin_back(3,i,iset),
2523 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2524 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2528 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2529 & wfrag_back(3,i,iset),
2530 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2531 if(me.eq.king.or..not.out1file)
2532 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2533 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2539 C---------------------------------------------------------------------------
2540 subroutine read_afminp
2542 include 'DIMENSIONS'
2546 include 'COMMON.SETUP'
2547 include 'COMMON.CONTROL'
2548 include 'COMMON.CHAIN'
2549 include 'COMMON.IOUNITS'
2550 include 'COMMON.SBRIDGE'
2551 character*320 afmcard
2553 c print *, "wchodze"
2554 call card_concat(afmcard)
2555 call readi(afmcard,"BEG",afmbeg,0)
2556 call readi(afmcard,"END",afmend,0)
2557 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2558 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2559 c print *,'FORCE=' ,forceAFMconst
2560 CCCC NOW PROPERTIES FOR AFM
2563 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2565 distafminit=dsqrt(distafminit)
2566 c print *,'initdist',distafminit
2569 c-------------------------------------------------------------------------------
2570 subroutine read_saxs_constr
2572 include 'DIMENSIONS'
2576 include 'COMMON.SETUP'
2577 include 'COMMON.CONTROL'
2578 include 'COMMON.SAXS'
2579 include 'COMMON.CHAIN'
2580 include 'COMMON.IOUNITS'
2581 include 'COMMON.SBRIDGE'
2582 double precision cm(3),cnorm
2585 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2587 if (saxs_mode.eq.0) then
2588 c SAXS distance distribution
2590 read(inp,*) distsaxs(i),Psaxs(i)
2594 Cnorm = Cnorm + Psaxs(i)
2596 write (iout,*) "Cnorm",Cnorm
2598 Psaxs(i)=Psaxs(i)/Cnorm
2600 write (iout,*) "Normalized distance distribution from SAXS"
2602 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2606 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2608 write (iout,*) "Wsaxs0",Wsaxs0
2612 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2619 cm(j)=cm(j)+Csaxs(j,i)
2627 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2630 write (iout,*) "SAXS sphere coordinates"
2632 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2638 c-------------------------------------------------------------------------------
2639 subroutine read_dist_constr
2641 include 'DIMENSIONS'
2645 include 'COMMON.SETUP'
2646 include 'COMMON.CONTROL'
2647 include 'COMMON.CHAIN'
2648 include 'COMMON.IOUNITS'
2649 include 'COMMON.SBRIDGE'
2650 include 'COMMON.INTERACT'
2651 integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
2652 integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
2653 double precision wfrag_(100),wpair_(1000)
2654 double precision ddjk,dist,dist_cut,fordepthmax
2655 character*5000 controlcard
2656 logical normalize,next
2658 double precision scal_bfac
2659 double precision xlink(4,0:4) /
2661 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2662 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2663 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2664 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2665 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2666 c print *, "WCHODZE"
2667 write (iout,*) "Calling read_dist_constr"
2668 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2670 restr_on_coord=.false.
2679 call card_concat(controlcard)
2680 next = index(controlcard,"NEXT").gt.0
2681 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2682 write (iout,*) "restr_type",restr_type
2683 call readi(controlcard,"NFRAG",nfrag_,0)
2684 call readi(controlcard,"NFRAG",nfrag_,0)
2685 call readi(controlcard,"NPAIR",npair_,0)
2686 call readi(controlcard,"NDIST",ndist_,0)
2687 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2688 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2689 if (restr_type.eq.10)
2690 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2691 if (restr_type.eq.12)
2692 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2693 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2694 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2695 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2696 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2697 normalize = index(controlcard,"NORMALIZE").gt.0
2698 write (iout,*) "WBOLTZD",wboltzd
2699 write (iout,*) "SCAL_PEAK",scal_peak
2700 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2701 write (iout,*) "IFRAG"
2703 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2705 write (iout,*) "IPAIR"
2707 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2709 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2711 & "Distance restraints as generated from reference structure"
2713 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2714 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2715 & ifrag_(2,i)=nstart_sup+nsup-1
2716 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2718 if (wfrag_(i).eq.0.0d0) cycle
2719 do j=ifrag_(1,i),ifrag_(2,i)-1
2720 do k=j+1,ifrag_(2,i)
2721 c write (iout,*) "j",j," k",k
2723 if (restr_type.eq.1) then
2729 forcon(nhpb)=wfrag_(i)
2730 else if (constr_dist.eq.2) then
2731 if (ddjk.le.dist_cut) then
2737 forcon(nhpb)=wfrag_(i)
2739 else if (restr_type.eq.3) then
2745 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2748 if (.not.out1file .or. me.eq.king)
2749 & write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
2750 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2752 write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
2753 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2759 if (wpair_(i).eq.0.0d0) cycle
2767 do j=ifrag_(1,ii),ifrag_(2,ii)
2768 do k=ifrag_(1,jj),ifrag_(2,jj)
2770 if (restr_type.eq.1) then
2776 forcon(nhpb)=wpair_(i)
2777 else if (constr_dist.eq.2) then
2778 if (ddjk.le.dist_cut) then
2784 forcon(nhpb)=wpair_(i)
2786 else if (restr_type.eq.3) then
2792 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2795 if (.not.out1file .or. me.eq.king)
2796 & write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
2797 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2799 write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
2800 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2807 write (iout,*) "Distance restraints as read from input"
2809 if (restr_type.eq.12) then
2810 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2811 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2812 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2813 & fordepth_peak(nhpb_peak+1),npeak
2814 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2815 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2816 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2817 c & fordepth_peak(nhpb_peak+1),npeak
2818 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2819 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2820 nhpb_peak=nhpb_peak+1
2821 irestr_type_peak(nhpb_peak)=12
2822 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2825 if (.not.out1file .or. me.eq.king)
2826 & write (iout,'(a,5i6,2f8.2,2f10.5,i5)') "+dist.restr ",
2827 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2828 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2829 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2830 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2832 write (iout,'(a,5i6,2f8.2,2f10.5,i5)') "+dist.restr ",
2833 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2834 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2835 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2836 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2838 if (ibecarb_peak(nhpb_peak).eq.3) then
2839 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2840 else if (ibecarb_peak(nhpb_peak).eq.2) then
2841 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2842 else if (ibecarb_peak(nhpb_peak).eq.1) then
2843 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2844 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2846 else if (restr_type.eq.11) then
2847 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2848 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2849 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2850 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2852 irestr_type(nhpb)=11
2854 if (.not.out1file .or. me.eq.king)
2855 & write (iout,'(a,4i6,2f8.2,2f10.5,i5)') "+dist.restr ",
2856 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2857 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2859 write (iout,'(a,4i6,2f8.2,2f10.5,i5)') "+dist.restr ",
2860 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2861 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2863 c if (ibecarb(nhpb).gt.0) then
2864 c ihpb(nhpb)=ihpb(nhpb)+nres
2865 c jhpb(nhpb)=jhpb(nhpb)+nres
2867 if (ibecarb(nhpb).eq.3) then
2868 ihpb(nhpb)=ihpb(nhpb)+nres
2869 else if (ibecarb(nhpb).eq.2) then
2870 ihpb(nhpb)=ihpb(nhpb)+nres
2871 else if (ibecarb(nhpb).eq.1) then
2872 ihpb(nhpb)=ihpb(nhpb)+nres
2873 jhpb(nhpb)=jhpb(nhpb)+nres
2875 else if (restr_type.eq.10) then
2876 c Cross-lonk Markov-like potential
2877 call card_concat(controlcard)
2878 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2879 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2881 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2882 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2883 if (index(controlcard,"ZL").gt.0) then
2885 else if (index(controlcard,"ADH").gt.0) then
2887 else if (index(controlcard,"PDH").gt.0) then
2889 else if (index(controlcard,"DSS").gt.0) then
2894 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2895 & xlink(1,link_type))
2896 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2897 & xlink(2,link_type))
2898 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2899 & xlink(3,link_type))
2900 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2901 & xlink(4,link_type))
2902 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2903 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2904 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2905 if (forcon(nhpb+1).le.0.0d0 .or.
2906 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2908 irestr_type(nhpb)=10
2909 if (ibecarb(nhpb).eq.3) then
2910 jhpb(nhpb)=jhpb(nhpb)+nres
2911 else if (ibecarb(nhpb).eq.2) then
2912 ihpb(nhpb)=ihpb(nhpb)+nres
2913 else if (ibecarb(nhpb).eq.1) then
2914 ihpb(nhpb)=ihpb(nhpb)+nres
2915 jhpb(nhpb)=jhpb(nhpb)+nres
2918 if (.not.out1file .or. me.eq.king)
2919 & write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
2920 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2921 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2924 write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
2925 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2926 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2931 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2932 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2933 if (forcon(nhpb+1).gt.0.0d0) then
2935 if (dhpb1(nhpb).eq.0.0d0) then
2940 if (ibecarb(nhpb).eq.3) then
2941 jhpb(nhpb)=jhpb(nhpb)+nres
2942 else if (ibecarb(nhpb).eq.2) then
2943 ihpb(nhpb)=ihpb(nhpb)+nres
2944 else if (ibecarb(nhpb).eq.1) then
2945 ihpb(nhpb)=ihpb(nhpb)+nres
2946 jhpb(nhpb)=jhpb(nhpb)+nres
2948 if (dhpb(nhpb).eq.0.0d0)
2949 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2952 if (.not.out1file .or. me.eq.king)
2953 & write (iout,'(a,4i6,f8.2,f10.1)') "+dist.restr ",
2954 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2956 write (iout,'(a,4i6,f8.2,f10.1)') "+dist.restr ",
2957 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2960 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2961 C if (forcon(nhpb+1).gt.0.0d0) then
2963 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2966 if (restr_type.eq.4) then
2967 write (iout,*) "The BFAC array"
2969 write (iout,'(i5,f10.5)') i,bfac(i)
2972 if (itype(i).eq.ntyp1) cycle
2974 if (itype(j).eq.ntyp1) cycle
2975 if (itype(i).eq.10) then
2980 if (itype(j).eq.10) then
2990 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2993 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2994 ihpb(nhpb)=i+nres*ii
2995 jhpb(nhpb)=j+nres*jj
2996 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2998 if (.not.out1file .or. me.eq.king) then
2999 write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
3000 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
3001 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
3005 write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
3006 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
3007 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
3017 if (restr_type.eq.5) then
3018 restr_on_coord=.true.
3020 if (itype(i).eq.ntyp1) cycle
3021 bfac(i)=(scal_bfac/bfac(i))**2
3030 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
3031 & fordepthmax=fordepth(i)
3034 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
3039 c-------------------------------------------------------------------------------
3040 subroutine read_constr_homology
3042 include 'DIMENSIONS'
3046 include 'COMMON.SETUP'
3047 include 'COMMON.CONTROL'
3048 include 'COMMON.HOMOLOGY'
3049 include 'COMMON.CHAIN'
3050 include 'COMMON.IOUNITS'
3052 include 'COMMON.QRESTR'
3053 include 'COMMON.GEO'
3054 include 'COMMON.INTERACT'
3055 include 'COMMON.NAMES'
3057 c For new homol impl
3059 include 'COMMON.VAR'
3062 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
3064 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
3065 c & sigma_odl_temp(maxres,maxres,max_template)
3067 character*24 model_ki_dist, model_ki_angle
3068 character*500 controlcard
3069 integer ki,i,ii,j,k,l,ii_in_use(maxdim_cont),i_tmp,idomain_tmp,
3070 & irec,ik,iistart,nres_temp
3073 logical liiflag,lfirst
3076 c FP - Nov. 2014 Temporary specifications for new vars
3078 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
3080 double precision, dimension (max_template,maxres) :: rescore
3081 double precision, dimension (max_template,maxres) :: rescore2
3082 double precision, dimension (max_template,maxres) :: rescore3
3083 double precision distal
3084 character*24 tpl_k_rescore
3085 character*256 pdbfile
3086 c -----------------------------------------------------------------
3087 c Reading multiple PDB ref structures and calculation of retraints
3088 c not using pre-computed ones stored in files model_ki_{dist,angle}
3090 c -----------------------------------------------------------------
3093 c Alternative: reading from input
3094 call card_concat(controlcard)
3095 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
3096 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
3097 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
3098 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
3099 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3100 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3101 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3102 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3103 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3104 if(.not.read2sigma.and.start_from_model) then
3105 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3106 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3107 start_from_model=.false.
3108 iranconf=(indpdb.le.0)
3110 if(start_from_model .and. (me.eq.king .or. .not. out1file))
3111 & write(iout,*) 'START_FROM_MODELS is ON'
3112 c if(start_from_model .and. rest) then
3113 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3114 c write(iout,*) 'START_FROM_MODELS is OFF'
3115 c write(iout,*) 'remove restart keyword from input'
3118 if (start_from_model) nmodel_start=constr_homology
3119 if (homol_nset.gt.1)then
3120 call card_concat(controlcard)
3121 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3122 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3123 c write(iout,*) "iset homology_weight "
3125 write(iout,*) i,waga_homology(i)
3128 iset=mod(kolor,homol_nset)+1
3131 waga_homology(1)=1.0
3134 cd write (iout,*) "nnt",nnt," nct",nct
3137 if (read_homol_frag) then
3138 call read_klapaucjusz
3144 c write(iout,*) 'nnt=',nnt,'nct=',nct
3147 c do k=1,constr_homology
3161 do k=1,constr_homology
3163 read(inp,'(a)') pdbfile
3164 pdbfiles_chomo(k)=pdbfile
3165 if(me.eq.king .or. .not. out1file)
3166 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3167 & pdbfile(:ilen(pdbfile))
3168 open(ipdbin,file=pdbfile,status='old',err=33)
3170 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3171 & pdbfile(:ilen(pdbfile))
3174 c print *,'Begin reading pdb data'
3176 c Files containing res sim or local scores (former containing sigmas)
3179 write(kic2,'(bz,i2.2)') k
3181 tpl_k_rescore="template"//kic2//".sco"
3185 if (read2sigma) then
3186 call readpdb_template(k)
3193 c Distance restraints
3196 C Copy the coordinates from reference coordinates (?)
3197 do i=1,2*nres_chomo(k)
3200 c write (iout,*) "c(",j,i,") =",c(j,i)
3204 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3206 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3207 open (ientin,file=tpl_k_rescore,status='old')
3208 if (nnt.gt.1) rescore(k,1)=0.0d0
3209 do irec=nnt,nct ! loop for reading res sim
3210 if (read2sigma) then
3211 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3212 & rescore3_tmp,idomain_tmp
3214 idomain(k,i_tmp)=idomain_tmp
3215 rescore(k,i_tmp)=rescore_tmp
3216 rescore2(k,i_tmp)=rescore2_tmp
3217 rescore3(k,i_tmp)=rescore3_tmp
3218 if (.not. out1file .or. me.eq.king)
3219 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3220 & i_tmp,rescore2_tmp,rescore_tmp,
3221 & rescore3_tmp,idomain_tmp
3224 read (ientin,*,end=1401) rescore_tmp
3226 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3227 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3228 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3233 if (waga_dist.ne.0.0d0) then
3241 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3242 c write (iout,*) k,i,j,distal,dist2_cut
3244 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3245 & .and. distal.le.dist2_cut ) then
3251 c write (iout,*) "k",k
3252 c write (iout,*) "i",i," j",j," constr_homology",
3257 if (read2sigma) then
3260 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3262 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3263 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3264 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3266 if (odl(k,ii).le.dist_cut) then
3267 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3270 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3271 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3273 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3274 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3278 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3281 c l_homo(k,ii)=.false.
3287 c write (iout,*) "Distance restraints set"
3290 c Theta, dihedral and SC retraints
3292 if (waga_angle.gt.0.0d0) then
3293 c open (ientin,file=tpl_k_sigma_dih,status='old')
3294 c do irec=1,maxres-3 ! loop for reading sigma_dih
3295 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3296 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3297 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3298 c & sigma_dih(k,i+nnt-1)
3303 if (idomain(k,i).eq.0) then
3307 dih(k,i)=phiref(i) ! right?
3308 c read (ientin,*) sigma_dih(k,i) ! original variant
3309 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3310 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3311 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3312 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3313 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3315 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3316 & rescore(k,i-2)+rescore(k,i-3))/4.0
3317 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3318 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3319 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3320 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3321 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3322 c Instead of res sim other local measure of b/b str reliability possible
3323 if (sigma_dih(k,i).ne.0)
3324 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3325 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3329 c write (iout,*) "Dihedral angle restraints set"
3332 if (waga_theta.gt.0.0d0) then
3333 c open (ientin,file=tpl_k_sigma_theta,status='old')
3334 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3335 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3336 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3337 c & sigma_theta(k,i+nnt-1)
3342 do i = nnt+2,nct ! right? without parallel.
3343 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3344 c do i=ithet_start,ithet_end ! with FG parallel.
3345 if (idomain(k,i).eq.0) then
3346 sigma_theta(k,i)=0.0
3349 thetatpl(k,i)=thetaref(i)
3350 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3351 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3352 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3353 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3354 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3355 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3356 & rescore(k,i-2))/3.0
3357 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3358 if (sigma_theta(k,i).ne.0)
3359 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3361 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3362 c rescore(k,i-2) ! right expression ?
3363 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3366 c write (iout,*) "Angle restraints set"
3369 if (waga_d.gt.0.0d0) then
3370 c open (ientin,file=tpl_k_sigma_d,status='old')
3371 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3372 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3373 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3374 c & sigma_d(k,i+nnt-1)
3378 do i = nnt,nct ! right? without parallel.
3379 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3380 c do i=loc_start,loc_end ! with FG parallel.
3381 if (itype(i).eq.10) cycle
3382 if (idomain(k,i).eq.0 ) then
3389 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3390 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3391 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3392 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3393 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3394 if (sigma_d(k,i).ne.0)
3395 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3397 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3398 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3399 c read (ientin,*) sigma_d(k,i) ! 1st variant
3403 c write (iout,*) "SC restraints set"
3406 c remove distance restraints not used in any model from the list
3407 c shift data in all arrays
3409 c write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
3410 if (waga_dist.ne.0.0d0) then
3417 c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3418 c & .and. distal.le.dist2_cut ) then
3419 c write (iout,*) "i",i," j",j," ii",ii
3421 if (ii_in_use(ii).eq.0.and.liiflag.or.
3422 & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3429 if(i10.eq.lim_odl) i10=i10+1
3431 ires_homo(iistart+ki)=ires_homo(ki+i01)
3432 jres_homo(iistart+ki)=jres_homo(ki+i01)
3433 ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3434 do k=1,constr_homology
3435 odl(k,iistart+ki)=odl(k,ki+i01)
3436 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3437 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3440 iistart=iistart+i10-i01
3443 if (ii_in_use(ii).ne.0.and..not.liiflag) then
3451 c write (iout,*) "Removing distances completed"
3453 endif ! .not. klapaucjusz
3455 if (constr_homology.gt.0) call homology_partition
3456 c write (iout,*) "After homology_partition"
3458 if (constr_homology.gt.0) call init_int_table
3459 c write (iout,*) "After init_int_table"
3461 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3462 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3466 if (.not.out_template_restr) return
3467 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3468 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3469 write (iout,*) "Distance restraints from templates"
3471 write(iout,'(3i7,100(2f8.2,1x,l1,4x))')
3472 & ii,ires_homo(ii),jres_homo(ii),
3473 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3474 & ki=1,constr_homology)
3476 write (iout,*) "Dihedral angle restraints from templates"
3478 write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3479 & (rad2deg*dih(ki,i),
3480 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3482 write (iout,*) "Virtual-bond angle restraints from templates"
3484 write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3485 & (rad2deg*thetatpl(ki,i),
3486 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3488 write (iout,*) "SC restraints from templates"
3490 write(iout,'(i7,100(4f8.2,4x))') i,
3491 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3492 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3495 c -----------------------------------------------------------------
3498 c----------------------------------------------------------------------
3500 subroutine flush(iu)
3505 subroutine flush(iu)
3510 c------------------------------------------------------------------------------
3511 subroutine copy_to_tmp(source)
3513 include "DIMENSIONS"
3514 include "COMMON.IOUNITS"
3515 character*(*) source
3516 character* 256 tmpfile
3520 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3521 inquire(file=tmpfile,exist=ex)
3523 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3524 & " to temporary directory..."
3525 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3526 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3530 c------------------------------------------------------------------------------
3531 subroutine move_from_tmp(source)
3533 include "DIMENSIONS"
3534 include "COMMON.IOUNITS"
3535 character*(*) source
3538 write (*,*) "Moving ",source(:ilen(source)),
3539 & " from temporary directory to working directory"
3540 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3541 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3544 c------------------------------------------------------------------------------
3545 subroutine random_init(seed)
3547 C Initialize random number generator
3550 include 'DIMENSIONS'
3553 logical OKRandom, prng_restart
3555 integer iseed_array(4)
3556 integer error_msg,ierr
3558 include 'COMMON.IOUNITS'
3559 include 'COMMON.TIME1'
3560 include 'COMMON.THREAD'
3561 include 'COMMON.SBRIDGE'
3562 include 'COMMON.CONTROL'
3563 include 'COMMON.MCM'
3564 include 'COMMON.MAP'
3565 include 'COMMON.HEADER'
3566 include 'COMMON.CSA'
3567 include 'COMMON.CHAIN'
3568 include 'COMMON.MUCA'
3570 include 'COMMON.FFIELD'
3571 include 'COMMON.SETUP'
3573 double precision seed,ran_number
3574 iseed=-dint(dabs(seed))
3575 if (iseed.eq.0) then
3576 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3577 & 'Random seed undefined. The program will stop.'
3578 write (*,'(/80(1h*)/20x,a/80(1h*))')
3579 & 'Random seed undefined. The program will stop.'
3581 call mpi_finalize(mpi_comm_world,ierr)
3583 stop 'Bad random seed.'
3586 if (fg_rank.eq.0) then
3590 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3591 OKRandom = prng_restart(me,iseed)
3594 tmp=65536.0d0**(4-i)
3595 iseed_array(i) = dint(seed/tmp)
3596 seed=seed-iseed_array(i)*tmp
3599 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3600 & (iseed_array(i),i=1,4)
3601 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3602 & (iseed_array(i),i=1,4)
3603 OKRandom = prng_restart(me,iseed_array)
3606 c r1 = prng_next(me)
3607 r1=ran_number(0.0D0,1.0D0)
3609 & write (iout,*) 'ran_num',r1
3610 if (r1.lt.0.0d0) OKRandom=.false.
3612 if (.not.OKRandom) then
3613 write (iout,*) 'PRNG IS NOT WORKING!!!'
3614 print *,'PRNG IS NOT WORKING!!!'
3617 call mpi_abort(mpi_comm_world,error_msg,ierr)
3620 write (iout,*) 'too many processors for parallel prng'
3621 write (*,*) 'too many processors for parallel prng'
3629 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3633 c----------------------------------------------------------------------
3634 subroutine read_klapaucjusz
3636 include 'DIMENSIONS'
3640 include 'COMMON.SETUP'
3641 include 'COMMON.CONTROL'
3642 include 'COMMON.HOMOLOGY'
3643 include 'COMMON.CHAIN'
3644 include 'COMMON.IOUNITS'
3646 include 'COMMON.GEO'
3647 include 'COMMON.INTERACT'
3648 include 'COMMON.NAMES'
3649 character*256 fragfile
3650 integer ninclust(maxclust),inclust(max_template,maxclust),
3651 & nresclust(maxclust),iresclust(maxres,maxclust),nclust
3654 character*24 model_ki_dist, model_ki_angle
3655 character*500 controlcard
3656 integer ki, i, j, jj,k, l, ii_in_use(maxdim_cont),i_tmp,
3658 & ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,
3659 & i01,i10,nnt_chain,nct_chain
3660 integer itype_temp(maxres)
3661 double precision distal
3662 logical lprn /.true./
3666 logical liiflag,lfirst
3669 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3670 double precision, dimension (max_template,maxres) :: rescore
3671 double precision, dimension (max_template,maxres) :: rescore2
3672 character*24 tpl_k_rescore
3673 character*256 pdbfile
3676 c For new homol impl
3678 include 'COMMON.VAR'
3680 c write (iout,*) "READ_KLAPAUCJUSZ"
3681 c print *,"READ_KLAPAUCJUSZ"
3683 call getenv("FRAGFILE",fragfile)
3684 write (iout,*) "Opening", fragfile
3686 open(ientin,file=fragfile,status="old",err=10)
3687 c write (iout,*) " opened"
3700 c write (iout,*) "Entering loop"
3703 DO IGR = 1,NCHAIN_GROUP
3705 c write (iout,*) "igr",igr
3707 read(ientin,*) constr_homology,nclust
3709 if (start_from_model) then
3710 nmodel_start=constr_homology
3717 ichain=iequiv(1,igr)
3718 nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
3719 nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
3720 c write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
3723 do k=1,constr_homology
3724 read(ientin,'(a)') pdbfile
3725 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3726 & pdbfile(:ilen(pdbfile))
3727 pdbfiles_chomo(k)=pdbfile
3728 open(ipdbin,file=pdbfile,status='old',err=33)
3730 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3731 & pdbfile(:ilen(pdbfile))
3736 call readpdb_template(k)
3746 read(ientin,*) ninclust(i),nresclust(i)
3747 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3748 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3751 c Loop over clusters
3754 do ll = 1,ninclust(l)
3757 c write (iout,*) "l",l," ll",ll," k",k
3763 idomain(k,iresclust(i,l)+1) = 1
3765 idomain(k,iresclust(i,l)) = 1
3769 c Distance restraints
3772 C Copy the coordinates from reference coordinates (?)
3778 c write (iout,*) "c(",j,i,") =",c(j,i)
3781 call int_from_cart(.true.,.false.)
3782 call sc_loc_geom(.false.)
3784 thetaref(i)=theta(i)
3788 if (waga_dist.ne.0.0d0) then
3791 do i = nnt_chain,nct_chain-2
3798 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3799 c write (iout,*) k,i,j,distal,dist2_cut
3801 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3802 & .and. distal.le.dist2_cut ) then
3808 c write (iout,*) "k",k
3809 c write (iout,*) "i",i," j",j," constr_homology",
3811 ires_homo(ii)=i+chain_border1(1,igr)-1
3812 jres_homo(ii)=j+chain_border1(1,igr)-1
3814 if (read2sigma) then
3817 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3819 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3820 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3821 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3823 if (odl(k,ii).le.dist_cut) then
3824 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3827 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3828 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3830 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3831 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3835 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3838 c l_homo(k,ii)=.false.
3845 c Theta, dihedral and SC retraints
3847 if (waga_angle.gt.0.0d0) then
3848 do i = nnt_chain+3,nct_chain
3849 iii=i+chain_border1(1,igr)-1
3850 if (idomain(k,i).eq.0) then
3851 c sigma_dih(k,i)=0.0
3854 dih(k,iii)=phiref(i)
3856 & (rescore(k,i)+rescore(k,i-1)+
3857 & rescore(k,i-2)+rescore(k,i-3))/4.0
3858 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3859 c & " sigma_dihed",sigma_dih(k,i)
3860 if (sigma_dih(k,iii).ne.0)
3861 & sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
3866 if (waga_theta.gt.0.0d0) then
3867 do i = nnt_chain+2,nct_chain
3868 iii=i+chain_border1(1,igr)-1
3869 if (idomain(k,i).eq.0) then
3870 c sigma_theta(k,i)=0.0
3873 thetatpl(k,iii)=thetaref(i)
3874 sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+
3875 & rescore(k,i-2))/3.0
3876 if (sigma_theta(k,iii).ne.0)
3877 & sigma_theta(k,iii)=1.0d0/
3878 & (sigma_theta(k,iii)*sigma_theta(k,iii))
3882 if (waga_d.gt.0.0d0) then
3883 do i = nnt_chain,nct_chain
3884 iii=i+chain_border1(1,igr)-1
3885 if (itype(i).eq.10) cycle
3886 if (idomain(k,i).eq.0 ) then
3890 xxtpl(k,iii)=xxref(i)
3891 yytpl(k,iii)=yyref(i)
3892 zztpl(k,iii)=zzref(i)
3893 sigma_d(k,iii)=rescore(k,i)
3894 if (sigma_d(k,iii).ne.0)
3895 & sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
3896 c if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3902 c remove distance restraints not used in any model from the list
3903 c shift data in all arrays
3905 c write (iout,*) "ii_old",ii_old
3906 if (waga_dist.ne.0.0d0) then
3908 write (iout,*) "Distance restraints from templates"
3910 write(iout,'(4i5,100(2f8.2,1x,l1,4x))')
3911 & iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii),
3912 & (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii),
3913 & ki=1,constr_homology)
3919 do i=nnt_chain,nct_chain-2
3922 c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3923 c & .and. distal.le.dist2_cut ) then
3924 c write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
3926 if (ii_in_use(ii).eq.0.and.liiflag.or.
3927 & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3934 if(i10.eq.lim_odl) i10=i10+1
3936 ires_homo(iistart+ki)=ires_homo(ki+i01)
3937 jres_homo(iistart+ki)=jres_homo(ki+i01)
3938 ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3939 do k=1,constr_homology
3940 odl(k,iistart+ki)=odl(k,ki+i01)
3941 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3942 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3945 iistart=iistart+i10-i01
3948 if (ii_in_use(ii).ne.0.and..not.liiflag) then
3961 ichain=iequiv(i,igr)
3963 do j=nnt_chain,nct_chain
3964 jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
3965 do k=1,constr_homology
3967 sigma_dih(k,jj)=sigma_dih(k,j)
3968 thetatpl(k,jj)=thetatpl(k,j)
3969 sigma_theta(k,jj)=sigma_theta(k,j)
3970 xxtpl(k,jj)=xxtpl(k,j)
3971 yytpl(k,jj)=yytpl(k,j)
3972 zztpl(k,jj)=zztpl(k,j)
3973 sigma_d(k,jj)=sigma_d(k,j)
3977 jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr))
3978 c write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
3979 do j=ii_old+1,lim_odl
3980 ires_homo(j+lll)=ires_homo(j)+jj
3981 jres_homo(j+lll)=jres_homo(j)+jj
3982 do k=1,constr_homology
3983 odl(k,j+lll)=odl(k,j)
3984 sigma_odl(k,j+lll)=sigma_odl(k,j)
3985 l_homo(k,j+lll)=l_homo(k,j)
3996 if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2
4001 10 stop "Error in fragment file"