2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 if (modecalc.eq.8) then
33 inquire (file="fort.40",exist=file_exist)
34 if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
75 include 'COMMON.HEADER'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.INTERACT'
82 include 'COMMON.SETUP'
83 include 'COMMON.SPLITELE'
84 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
85 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
87 character*320 controlcard
92 read (INP,'(a)') titel
93 call card_concat(controlcard)
94 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
95 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
96 call reada(controlcard,'SEED',seed,0.0D0)
97 call random_init(seed)
98 C Set up the time limit (caution! The time must be input in minutes!)
99 read_cart=index(controlcard,'READ_CART').gt.0
100 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
101 write (iout,*) "constr_dist",constr_dist
102 call readi(controlcard,'NSAXS',nsaxs,0)
103 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
104 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
105 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
106 & SAXS_MODE," SCAL_RAD",scal_rad
107 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
108 call readi(controlcard,'SYM',symetr,1)
109 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
110 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
111 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
112 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
113 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
114 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
115 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
116 call reada(controlcard,'DRMS',drms,0.1D0)
117 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
118 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
119 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
120 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
121 write (iout,'(a,f10.1)')'DRMS = ',drms
122 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
123 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
125 call readi(controlcard,'NZ_START',nz_start,0)
126 call readi(controlcard,'NZ_END',nz_end,0)
127 call readi(controlcard,'IZ_SC',iz_sc,0)
129 safety = 60.0d0*safety
132 call reada(controlcard,"T_BATH",t_bath,300.0d0)
133 minim=(index(controlcard,'MINIMIZE').gt.0)
134 dccart=(index(controlcard,'CART').gt.0)
135 overlapsc=(index(controlcard,'OVERLAP').gt.0)
136 overlapsc=.not.overlapsc
137 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
138 searchsc=.not.searchsc
139 sideadd=(index(controlcard,'SIDEADD').gt.0)
140 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
141 outpdb=(index(controlcard,'PDBOUT').gt.0)
142 outmol2=(index(controlcard,'MOL2OUT').gt.0)
143 pdbref=(index(controlcard,'PDBREF').gt.0)
144 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
145 indpdb=index(controlcard,'PDBSTART')
146 extconf=(index(controlcard,'EXTCONF').gt.0)
147 AFMlog=(index(controlcard,'AFM'))
148 selfguide=(index(controlcard,'SELFGUIDE'))
149 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
150 call readi(controlcard,'IPRINT',iprint,0)
151 call readi(controlcard,'MAXGEN',maxgen,10000)
152 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
153 call readi(controlcard,"KDIAG",kdiag,0)
154 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
155 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
156 & write (iout,*) "RESCALE_MODE",rescale_mode
157 split_ene=index(controlcard,'SPLIT_ENE').gt.0
158 if (index(controlcard,'REGULAR').gt.0.0D0) then
159 call reada(controlcard,'WEIDIS',weidis,0.1D0)
163 if (index(controlcard,'CHECKGRAD').gt.0) then
165 if (index(controlcard,'CART').gt.0) then
167 elseif (index(controlcard,'CARINT').gt.0) then
172 elseif (index(controlcard,'THREAD').gt.0) then
174 call readi(controlcard,'THREAD',nthread,0)
175 if (nthread.gt.0) then
176 call reada(controlcard,'WEIDIS',weidis,0.1D0)
179 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
180 stop 'Error termination in Read_Control.'
182 else if (index(controlcard,'MCMA').gt.0) then
184 else if (index(controlcard,'MCEE').gt.0) then
186 else if (index(controlcard,'MULTCONF').gt.0) then
188 else if (index(controlcard,'MAP').gt.0) then
190 call readi(controlcard,'MAP',nmap,0)
191 else if (index(controlcard,'CSA').gt.0) then
193 crc else if (index(controlcard,'ZSCORE').gt.0) then
195 crc ZSCORE is rm from UNRES, modecalc=9 is available
198 cfcm else if (index(controlcard,'MCMF').gt.0) then
200 else if (index(controlcard,'SOFTREG').gt.0) then
202 else if (index(controlcard,'CHECK_BOND').gt.0) then
204 else if (index(controlcard,'TEST').gt.0) then
206 else if (index(controlcard,'MD').gt.0) then
208 else if (index(controlcard,'RE ').gt.0) then
212 lmuca=index(controlcard,'MUCA').gt.0
213 call readi(controlcard,'MUCADYN',mucadyn,0)
214 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
215 if (lmuca .and. (me.eq.king .or. .not.out1file ))
217 write (iout,*) 'MUCADYN=',mucadyn
218 write (iout,*) 'MUCASMOOTH=',muca_smooth
221 iscode=index(controlcard,'ONE_LETTER')
222 indphi=index(controlcard,'PHI')
223 indback=index(controlcard,'BACK')
224 iranconf=index(controlcard,'RAND_CONF')
225 i2ndstr=index(controlcard,'USE_SEC_PRED')
226 gradout=index(controlcard,'GRADOUT').gt.0
227 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
228 C DISTCHAINMAX become obsolete for periodic boundry condition
229 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
230 C Reading the dimensions of box in x,y,z coordinates
231 call reada(controlcard,'BOXX',boxxsize,100.0d0)
232 call reada(controlcard,'BOXY',boxysize,100.0d0)
233 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
234 c Cutoff range for interactions
235 call reada(controlcard,"R_CUT",r_cut,15.0d0)
236 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
237 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
238 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
239 if (lipthick.gt.0.0d0) then
240 bordliptop=(boxzsize+lipthick)/2.0
241 bordlipbot=bordliptop-lipthick
243 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
244 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
245 buflipbot=bordlipbot+lipbufthick
246 bufliptop=bordliptop-lipbufthick
247 if ((lipbufthick*2.0d0).gt.lipthick)
248 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
250 c write(iout,*) "bordliptop=",bordliptop
251 c write(iout,*) "bordlipbot=",bordlipbot
252 c write(iout,*) "bufliptop=",bufliptop
253 c write(iout,*) "buflipbot=",buflipbot
256 if (me.eq.king .or. .not.out1file )
257 & write (iout,*) "DISTCHAINMAX",distchainmax
259 if(me.eq.king.or..not.out1file)
260 & write (iout,'(2a)') diagmeth(kdiag),
261 & ' routine used to diagonalize matrices.'
264 c--------------------------------------------------------------------------
265 subroutine read_REMDpar
269 implicit real*8 (a-h,o-z)
271 include 'COMMON.IOUNITS'
272 include 'COMMON.TIME1'
275 include 'COMMON.LANGEVIN'
277 include 'COMMON.LANGEVIN.lang0'
279 include 'COMMON.INTERACT'
280 include 'COMMON.NAMES'
282 include 'COMMON.REMD'
283 include 'COMMON.CONTROL'
284 include 'COMMON.SETUP'
286 character*320 controlcard
287 character*3200 controlcard1
288 integer iremd_m_total
290 if(me.eq.king.or..not.out1file)
291 & write (iout,*) "REMD setup"
293 call card_concat(controlcard)
294 call readi(controlcard,"NREP",nrep,3)
295 call readi(controlcard,"NSTEX",nstex,1000)
296 call reada(controlcard,"RETMIN",retmin,10.0d0)
297 call reada(controlcard,"RETMAX",retmax,1000.0d0)
298 mremdsync=(index(controlcard,'SYNC').gt.0)
299 call readi(controlcard,"NSYN",i_sync_step,100)
300 restart1file=(index(controlcard,'REST1FILE').gt.0)
301 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
302 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
303 if(max_cache_traj_use.gt.max_cache_traj)
304 & max_cache_traj_use=max_cache_traj
305 if(me.eq.king.or..not.out1file) then
306 cd if (traj1file) then
307 crc caching is in testing - NTWX is not ignored
308 cd write (iout,*) "NTWX value is ignored"
309 cd write (iout,*) " trajectory is stored to one file by master"
310 cd write (iout,*) " before exchange at NSTEX intervals"
312 write (iout,*) "NREP= ",nrep
313 write (iout,*) "NSTEX= ",nstex
314 write (iout,*) "SYNC= ",mremdsync
315 write (iout,*) "NSYN= ",i_sync_step
316 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
319 if (index(controlcard,'TLIST').gt.0) then
321 call card_concat(controlcard1)
322 read(controlcard1,*) (remd_t(i),i=1,nrep)
323 if(me.eq.king.or..not.out1file)
324 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
327 if (index(controlcard,'MLIST').gt.0) then
329 call card_concat(controlcard1)
330 read(controlcard1,*) (remd_m(i),i=1,nrep)
331 if(me.eq.king.or..not.out1file) then
332 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
335 iremd_m_total=iremd_m_total+remd_m(i)
337 write (iout,*) 'Total number of replicas ',iremd_m_total
340 if(me.eq.king.or..not.out1file)
341 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
344 c--------------------------------------------------------------------------
345 subroutine read_MDpar
349 implicit real*8 (a-h,o-z)
351 include 'COMMON.IOUNITS'
352 include 'COMMON.TIME1'
355 include 'COMMON.LANGEVIN'
357 include 'COMMON.LANGEVIN.lang0'
359 include 'COMMON.INTERACT'
360 include 'COMMON.NAMES'
362 include 'COMMON.SETUP'
363 include 'COMMON.CONTROL'
364 include 'COMMON.SPLITELE'
366 character*320 controlcard
368 call card_concat(controlcard)
369 call readi(controlcard,"NSTEP",n_timestep,1000000)
370 call readi(controlcard,"NTWE",ntwe,100)
371 call readi(controlcard,"NTWX",ntwx,1000)
372 call reada(controlcard,"DT",d_time,1.0d-1)
373 call reada(controlcard,"DVMAX",dvmax,2.0d1)
374 call reada(controlcard,"DAMAX",damax,1.0d1)
375 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
376 call readi(controlcard,"LANG",lang,0)
377 RESPA = index(controlcard,"RESPA") .gt. 0
378 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
379 ntime_split0=ntime_split
380 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
381 ntime_split0=ntime_split
382 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
383 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
384 rest = index(controlcard,"REST").gt.0
385 tbf = index(controlcard,"TBF").gt.0
386 usampl = index(controlcard,"USAMPL").gt.0
388 mdpdb = index(controlcard,"MDPDB").gt.0
389 call reada(controlcard,"T_BATH",t_bath,300.0d0)
390 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
391 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
392 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
393 if (count_reset_moment.eq.0) count_reset_moment=1000000000
394 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
395 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
396 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
397 if (count_reset_vel.eq.0) count_reset_vel=1000000000
398 large = index(controlcard,"LARGE").gt.0
399 print_compon = index(controlcard,"PRINT_COMPON").gt.0
400 rattle = index(controlcard,"RATTLE").gt.0
401 preminim = index(controlcard,"PREMINIM").gt.0
403 dccart=(index(controlcard,'CART').gt.0)
406 c if performing umbrella sampling, fragments constrained are read from the fragment file
412 if(me.eq.king.or..not.out1file) then
414 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
416 write (iout,'(a)') "The units are:"
417 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
418 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
419 & " acceleration: angstrom/(48.9 fs)**2"
420 write (iout,'(a)') "energy: kcal/mol, temperature: K"
422 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
423 write (iout,'(a60,f10.5,a)')
424 & "Initial time step of numerical integration:",d_time,
426 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
428 write (iout,'(2a,i4,a)')
429 & "A-MTS algorithm used; initial time step for fast-varying",
430 & " short-range forces split into",ntime_split," steps."
431 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
432 & r_cut," lambda",rlamb
434 write (iout,'(2a,f10.5)')
435 & "Maximum acceleration threshold to reduce the time step",
436 & "/increase split number:",damax
437 write (iout,'(2a,f10.5)')
438 & "Maximum predicted energy drift to reduce the timestep",
439 & "/increase split number:",edriftmax
440 write (iout,'(a60,f10.5)')
441 & "Maximum velocity threshold to reduce velocities:",dvmax
442 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
443 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
444 if (rattle) write (iout,'(a60)')
445 & "Rattle algorithm used to constrain the virtual bonds"
446 if (preminim .or. iranconf.gt.0) then
448 & "Initial structure will be energy-minimized"
453 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
454 call reada(controlcard,"RWAT",rwat,1.4d0)
455 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
456 surfarea=index(controlcard,"SURFAREA").gt.0
457 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
458 if(me.eq.king.or..not.out1file)then
459 write (iout,'(/a,$)') "Langevin dynamics calculation"
462 & " with direct integration of Langevin equations"
463 else if (lang.eq.2) then
464 write (iout,'(a/)') " with TINKER stochasic MD integrator"
465 else if (lang.eq.3) then
466 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
467 else if (lang.eq.4) then
468 write (iout,'(a/)') " in overdamped mode"
470 write (iout,'(//a,i5)')
471 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
474 write (iout,'(a60,f10.5)') "Temperature:",t_bath
475 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
476 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
477 write (iout,'(a60,f10.5)')
478 & "Scaling factor of the friction forces:",scal_fric
479 if (surfarea) write (iout,'(2a,i10,a)')
480 & "Friction coefficients will be scaled by solvent-accessible",
481 & " surface area every",reset_fricmat," steps."
483 c Calculate friction coefficients and bounds of stochastic forces
484 eta=6*pi*cPoise*etawat
485 if(me.eq.king.or..not.out1file)
486 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
488 gamp=scal_fric*(pstok+rwat)*eta
489 stdfp=dsqrt(2*Rb*t_bath/d_time)
491 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
492 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
494 if(me.eq.king.or..not.out1file)then
495 write (iout,'(/2a/)')
496 & "Radii of site types and friction coefficients and std's of",
497 & " stochastic forces of fully exposed sites"
498 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
500 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
501 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
505 if(me.eq.king.or..not.out1file)then
506 write (iout,'(a)') "Berendsen bath calculation"
507 write (iout,'(a60,f10.5)') "Temperature:",t_bath
508 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
510 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
511 & count_reset_moment," steps"
513 & write (iout,'(a,i10,a)')
514 & "Velocities will be reset at random every",count_reset_vel,
518 if(me.eq.king.or..not.out1file)
519 & write (iout,'(a31)') "Microcanonical mode calculation"
521 if(me.eq.king.or..not.out1file)then
522 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
524 write(iout,*) "MD running with constraints."
525 write(iout,*) "Equilibration time ", eq_time, " mtus."
526 write(iout,*) "Constraining ", nfrag," fragments."
527 write(iout,*) "Length of each fragment, weight and q0:"
529 write (iout,*) "Set of restraints #",iset
531 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
532 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
534 write(iout,*) "constraints between ", npair, "fragments."
535 write(iout,*) "constraint pairs, weights and q0:"
537 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
538 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
540 write(iout,*) "angle constraints within ", nfrag_back,
541 & "backbone fragments."
542 write(iout,*) "fragment, weights:"
544 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
545 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
546 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
549 iset=mod(kolor,nset)+1
552 if(me.eq.king.or..not.out1file)
553 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
556 c------------------------------------------------------------------------------
559 C Read molecular data.
561 implicit real*8 (a-h,o-z)
567 include 'COMMON.IOUNITS'
570 include 'COMMON.INTERACT'
571 include 'COMMON.LOCAL'
572 include 'COMMON.NAMES'
573 include 'COMMON.CHAIN'
574 include 'COMMON.FFIELD'
575 include 'COMMON.SBRIDGE'
576 include 'COMMON.HEADER'
577 include 'COMMON.CONTROL'
578 include 'COMMON.DBASE'
579 include 'COMMON.THREAD'
580 include 'COMMON.CONTACTS'
581 include 'COMMON.TORCNSTR'
582 include 'COMMON.TIME1'
583 include 'COMMON.BOUNDS'
585 include 'COMMON.SETUP'
586 character*4 sequence(maxres)
588 double precision x(maxvar)
589 character*256 pdbfile
590 character*320 weightcard
591 character*80 weightcard_t,ucase
592 dimension itype_pdb(maxres)
593 common /pizda/ itype_pdb
594 logical seq_comp,fail
595 double precision energia(0:n_ene)
601 C Read weights of the subsequent energy terms.
602 call card_concat(weightcard)
603 call reada(weightcard,'WLONG',wlong,1.0D0)
604 call reada(weightcard,'WSC',wsc,wlong)
605 call reada(weightcard,'WSCP',wscp,wlong)
606 call reada(weightcard,'WELEC',welec,1.0D0)
607 call reada(weightcard,'WVDWPP',wvdwpp,welec)
608 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
609 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
610 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
611 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
612 call reada(weightcard,'WTURN3',wturn3,1.0D0)
613 call reada(weightcard,'WTURN4',wturn4,1.0D0)
614 call reada(weightcard,'WTURN6',wturn6,1.0D0)
615 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
616 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
617 call reada(weightcard,'WBOND',wbond,1.0D0)
618 call reada(weightcard,'WTOR',wtor,1.0D0)
619 call reada(weightcard,'WTORD',wtor_d,1.0D0)
620 call reada(weightcard,'WANG',wang,1.0D0)
621 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
622 call reada(weightcard,'SCAL14',scal14,0.4D0)
623 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
624 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
625 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
626 call reada(weightcard,'TEMP0',temp0,300.0d0)
627 call reada(weightcard,'WLT',wliptran,0.0D0)
628 call reada(weightcard,'WSAXS',wsaxs,1.0D0)
629 if (index(weightcard,'SOFT').gt.0) ipot=6
630 C 12/1/95 Added weight for the multi-body term WCORR
631 call reada(weightcard,'WCORRH',wcorr,1.0D0)
632 if (wcorr4.gt.0.0d0) wcorr=wcorr4
653 if(me.eq.king.or..not.out1file)
654 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
655 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
657 10 format (/'Energy-term weights (unscaled):'//
658 & 'WSCC= ',f10.6,' (SC-SC)'/
659 & 'WSCP= ',f10.6,' (SC-p)'/
660 & 'WELEC= ',f10.6,' (p-p electr)'/
661 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
662 & 'WBOND= ',f10.6,' (stretching)'/
663 & 'WANG= ',f10.6,' (bending)'/
664 & 'WSCLOC= ',f10.6,' (SC local)'/
665 & 'WTOR= ',f10.6,' (torsional)'/
666 & 'WTORD= ',f10.6,' (double torsional)'/
667 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
668 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
669 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
670 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
671 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
672 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
673 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
674 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
675 & 'WTURN6= ',f10.6,' (turns, 6th order)')
676 if(me.eq.king.or..not.out1file)then
677 if (wcorr4.gt.0.0d0) then
678 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
679 & 'between contact pairs of peptide groups'
680 write (iout,'(2(a,f5.3/))')
681 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
682 & 'Range of quenching the correlation terms:',2*delt_corr
683 else if (wcorr.gt.0.0d0) then
684 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
685 & 'between contact pairs of peptide groups'
687 write (iout,'(a,f8.3)')
688 & 'Scaling factor of 1,4 SC-p interactions:',scal14
689 write (iout,'(a,f8.3)')
690 & 'General scaling factor of SC-p interactions:',scalscp
692 r0_corr=cutoff_corr-delt_corr
694 aad(i,1)=scalscp*aad(i,1)
695 aad(i,2)=scalscp*aad(i,2)
696 bad(i,1)=scalscp*bad(i,1)
697 bad(i,2)=scalscp*bad(i,2)
699 call rescale_weights(t_bath)
700 if(me.eq.king.or..not.out1file)
701 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
702 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
704 22 format (/'Energy-term weights (scaled):'//
705 & 'WSCC= ',f10.6,' (SC-SC)'/
706 & 'WSCP= ',f10.6,' (SC-p)'/
707 & 'WELEC= ',f10.6,' (p-p electr)'/
708 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
709 & 'WBOND= ',f10.6,' (stretching)'/
710 & 'WANG= ',f10.6,' (bending)'/
711 & 'WSCLOC= ',f10.6,' (SC local)'/
712 & 'WTOR= ',f10.6,' (torsional)'/
713 & 'WTORD= ',f10.6,' (double torsional)'/
714 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
715 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
716 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
717 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
718 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
719 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
720 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
721 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
722 & 'WTURN6= ',f10.6,' (turns, 6th order)')
723 if(me.eq.king.or..not.out1file)
724 & write (iout,*) "Reference temperature for weights calculation:",
726 call reada(weightcard,"D0CM",d0cm,3.78d0)
727 call reada(weightcard,"AKCM",akcm,15.1d0)
728 call reada(weightcard,"AKTH",akth,11.0d0)
729 call reada(weightcard,"AKCT",akct,12.0d0)
730 call reada(weightcard,"V1SS",v1ss,-1.08d0)
731 call reada(weightcard,"V2SS",v2ss,7.61d0)
732 call reada(weightcard,"V3SS",v3ss,13.7d0)
733 call reada(weightcard,"EBR",ebr,-5.50D0)
734 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
736 dyn_ss_mask(i)=.false.
740 dyn_ssbond_ij(i,j)=1.0d300
743 call reada(weightcard,"HT",Ht,0.0D0)
745 ss_depth=ebr/wsc-0.25*eps(1,1)
746 Ht=Ht/wsc-0.25*eps(1,1)
747 akcm=akcm*wstrain/wsc
748 akth=akth*wstrain/wsc
749 akct=akct*wstrain/wsc
750 v1ss=v1ss*wstrain/wsc
751 v2ss=v2ss*wstrain/wsc
752 v3ss=v3ss*wstrain/wsc
754 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
757 if(me.eq.king.or..not.out1file) then
758 write (iout,*) "Parameters of the SS-bond potential:"
759 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
761 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
762 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
763 write (iout,*)" HT",Ht
764 print *,'indpdb=',indpdb,' pdbref=',pdbref
766 if (indpdb.gt.0 .or. pdbref) then
767 read(inp,'(a)') pdbfile
768 if(me.eq.king.or..not.out1file)
769 & write (iout,'(2a)') 'PDB data will be read from file ',
770 & pdbfile(:ilen(pdbfile))
771 open(ipdbin,file=pdbfile,status='old',err=33)
773 33 write (iout,'(a)') 'Error opening PDB file.'
776 c print *,'Begin reading pdb data'
785 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
786 & (crefjlee(j,i+nres),j=1,3)
789 c print *,'Finished reading pdb data'
790 if(me.eq.king.or..not.out1file)
791 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
792 & ' nstart_sup=',nstart_sup
794 itype_pdb(i)=itype(i)
798 nct=nstart_sup+nsup-1
799 call contact(.false.,ncont_ref,icont_ref,co)
802 if(me.eq.king.or..not.out1file)
803 & write(iout,*)'Adding sidechains'
807 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
810 do while (fail.and.nsi.le.maxsi)
811 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
815 c Calculalte only the coordinates of the current sidechain; no need to rebuild
817 call locate_side_chain(i)
818 if(fail) write(iout,*)'Adding sidechain failed for res ',
819 & i,' after ',nsi,' trials'
822 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
826 if (indpdb.eq.0) then
827 C Read sequence if not taken from the pdb file.
829 c print *,'nres=',nres
830 if (iscode.gt.0) then
831 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
833 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
835 C Convert sequence to numeric code
837 itype(i)=rescode(i,sequence(i),iscode)
839 C Assign initial virtual bond lengths
845 vbld(i+nres)=dsc(iabs(itype(i)))
846 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
847 c write (iout,*) "i",i," itype",itype(i),
848 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
852 c print '(20i4)',(itype(i),i=1,nres)
855 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
857 if (itype(i).eq.ntyp1) then
861 else if (iabs(itype(i+1)).ne.20) then
863 else if (iabs(itype(i)).ne.20) then
870 if(me.eq.king.or..not.out1file)then
871 write (iout,*) "ITEL"
873 write (iout,*) i,itype(i),itel(i)
875 print *,'Call Read_Bridge.'
878 C 8/13/98 Set limits to generating the dihedral angles
883 read (inp,*) ndih_constr
884 write (iout,*) "ndish_constr",ndih_constr
885 if (ndih_constr.gt.0) then
887 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
888 if(me.eq.king.or..not.out1file)then
890 & 'There are',ndih_constr,' constraints on phi angles.'
892 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
896 phi0(i)=deg2rad*phi0(i)
897 drange(i)=deg2rad*drange(i)
899 if(me.eq.king.or..not.out1file)
900 & write (iout,*) 'FTORS',ftors
903 phibound(1,ii) = phi0(i)-drange(i)
904 phibound(2,ii) = phi0(i)+drange(i)
911 write (iout,'(a)') 'Boundaries in phi angle sampling:'
913 write (iout,'(a3,i5,2f10.1)')
914 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
920 cd print *,'NNT=',NNT,' NCT=',NCT
921 if (itype(1).eq.ntyp1) nnt=2
922 if (itype(nres).eq.ntyp1) nct=nct-1
924 if(me.eq.king.or..not.out1file)
925 & write (iout,'(a,i3)') 'nsup=',nsup
927 if (nsup.le.(nct-nnt+1)) then
928 do i=0,nct-nnt+1-nsup
929 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
935 & 'Error - sequences to be superposed do not match.'
938 do i=0,nsup-(nct-nnt+1)
939 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
941 nstart_sup=nstart_sup+i
947 & 'Error - sequences to be superposed do not match.'
950 if (nsup.eq.0) nsup=nct-nnt
951 if (nstart_sup.eq.0) nstart_sup=nnt
952 if (nstart_seq.eq.0) nstart_seq=nnt
953 if(me.eq.king.or..not.out1file)
954 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
955 & ' nstart_seq=',nstart_seq
957 c--- Zscore rms -------
958 if (nz_start.eq.0) nz_start=nnt
959 if (nz_end.eq.0 .and. nsup.gt.0) then
961 else if (nz_end.eq.0) then
964 if(me.eq.king.or..not.out1file)then
965 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
966 write (iout,*) 'IZ_SC=',iz_sc
968 c----------------------
971 if (.not.pdbref) then
972 call read_angles(inp,*38)
974 38 write (iout,'(a)') 'Error reading reference structure.'
976 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
977 stop 'Error reading reference structure'
979 39 call chainbuild_extconf
981 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
991 call contact(.true.,ncont_ref,icont_ref,co)
993 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
995 if (constr_dist.gt.0) call read_dist_constr
996 c write (iout,*) "After read_dist_constr nhpb",nhpb
997 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
998 if(me.eq.king.or..not.out1file)
999 & write (iout,*) 'Contact order:',co
1001 if(me.eq.king.or..not.out1file)
1002 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1005 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1007 if(me.eq.king.or..not.out1file)
1008 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1009 & icont_ref(1,i),' ',
1010 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1014 write (iout,*) "calling read_saxs_consrtr",nsaxs
1015 if (nsaxs.gt.0) call read_saxs_constr
1018 if (constr_homology.gt.0) then
1019 call read_constr_homology
1020 if (indpdb.gt.0 .or. pdbref) then
1023 c(j,i)=crefjlee(j,i)
1024 cref(j,i,1)=crefjlee(j,i)
1029 write (iout,*) "Array C"
1031 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1032 & (c(j,i+nres),j=1,3)
1034 write (iout,*) "Array Cref"
1036 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),
1037 & (cref(j,i+nres,1),j=1,3)
1040 call int_from_cart1(.false.)
1041 call sc_loc_geom(.false.)
1043 thetaref(i)=theta(i)
1048 dc(j,i)=c(j,i+1)-c(j,i)
1049 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1054 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1055 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1063 if (nhpb.gt.0) call hpb_partition
1064 c write (iout,*) "After read_dist_constr nhpb",nhpb
1066 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1067 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1068 & modecalc.ne.10) then
1069 C If input structure hasn't been supplied from the PDB file read or generate
1071 if (iranconf.eq.0 .and. .not. extconf) then
1072 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1073 & write (iout,'(a)') 'Initial geometry will be read in.'
1075 read(inp,'(8f10.5)',end=36,err=36)
1076 & ((c(l,k),l=1,3),k=1,nres),
1077 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1078 write (iout,*) "Exit READ_CART"
1079 write (iout,'(8f10.5)')
1080 & ((c(l,k),l=1,3),k=1,nres),
1081 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1082 call int_from_cart1(.true.)
1083 write (iout,*) "Finish INT_TO_CART"
1086 dc(j,i)=c(j,i+1)-c(j,i)
1087 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1091 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1093 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1094 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1100 write (iout,*) "Calling read_ang"
1101 call read_angles(inp,*36)
1102 write (iout,*) "Calling chainbuild"
1103 call chainbuild_extconf
1106 36 write (iout,'(a)') 'Error reading angle file.'
1108 call mpi_finalize( MPI_COMM_WORLD,IERR )
1110 stop 'Error reading angle file.'
1112 else if (extconf) then
1113 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1114 & write (iout,'(a)') 'Extended chain initial geometry.'
1116 theta(i)=90d0*deg2rad
1119 phi(i)=180d0*deg2rad
1122 alph(i)=110d0*deg2rad
1125 omeg(i)=-120d0*deg2rad
1126 if (itype(i).le.0) omeg(i)=-omeg(i)
1128 c from old chainbuild
1130 C Define the origin and orientation of the coordinate system and locate the
1131 C first three CA's and SC(2).
1135 * Build the alpha-carbon chain.
1138 call locate_next_res(i)
1141 C First and last SC must coincide with the corresponding CA.
1145 dc_norm(j,nres+1)=0.0D0
1146 dc(j,nres+nres)=0.0D0
1147 dc_norm(j,nres+nres)=0.0D0
1149 c(j,nres+nres)=c(j,nres)
1152 C Define the origin and orientation of the coordinate system and locate the
1153 C first three CA's and SC(2).
1157 * Build the alpha-carbon chain.
1160 call locate_next_res(i)
1163 C First and last SC must coincide with the corresponding CA.
1167 dc_norm(j,nres+1)=0.0D0
1168 dc(j,nres+nres)=0.0D0
1169 dc_norm(j,nres+nres)=0.0D0
1171 c(j,nres+nres)=c(j,nres)
1176 if(me.eq.king.or..not.out1file)
1177 & write (iout,'(a)') 'Random-generated initial geometry.'
1181 if (me.eq.king .or. fg_rank.eq.0 .and. (
1182 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1186 call gen_rand_conf(itmp,*30)
1188 30 write (iout,*) 'Failed to generate random conformation',
1189 & ', itrial=',itrial
1190 write (*,*) 'Processor:',me,
1191 & ' Failed to generate random conformation',
1201 write (iout,'(a,i3,a)') 'Processor:',me,
1202 & ' error in generating random conformation.'
1203 write (*,'(a,i3,a)') 'Processor:',me,
1204 & ' error in generating random conformation.'
1207 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1212 & ' error in generating random conformation.'
1217 elseif (modecalc.eq.4) then
1218 read (inp,'(a)') intinname
1219 open (intin,file=intinname,status='old',err=333)
1220 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1221 & write (iout,'(a)') 'intinname',intinname
1222 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1224 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1226 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1228 stop 'Error opening angle file.'
1232 C Generate distance constraints, if the PDB structure is to be regularized.
1233 if (nthread.gt.0) then
1234 call read_threadbase
1237 if (me.eq.king .or. .not. out1file)
1239 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1240 write (iout,'(/a,i3,a)')
1241 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1242 write (iout,'(20i4)') (iss(i),i=1,ns)
1244 write(iout,*)"Running with dynamic disulfide-bond formation"
1246 write (iout,'(/a/)') 'Pre-formed links are:'
1252 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1253 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1259 if (ns.gt.0.and.dyn_ss) then
1263 forcon(i-nss)=forcon(i)
1270 dyn_ss_mask(iss(i))=.true.
1273 if (i2ndstr.gt.0) call secstrp2dihc
1274 c call geom_to_var(nvar,x)
1275 c call etotal(energia(0))
1276 c call enerprint(energia(0))
1277 c call briefout(0,etot)
1279 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1280 cd write (iout,'(a)') 'Variable list:'
1281 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1283 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1284 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1285 & 'Processor',myrank,': end reading molecular data.'
1289 c--------------------------------------------------------------------------
1290 logical function seq_comp(itypea,itypeb,length)
1292 integer length,itypea(length),itypeb(length)
1295 if (itypea(i).ne.itypeb(i)) then
1303 c-----------------------------------------------------------------------------
1304 subroutine read_bridge
1305 C Read information about disulfide bridges.
1306 implicit real*8 (a-h,o-z)
1307 include 'DIMENSIONS'
1311 include 'COMMON.IOUNITS'
1312 include 'COMMON.GEO'
1313 include 'COMMON.VAR'
1314 include 'COMMON.INTERACT'
1315 include 'COMMON.LOCAL'
1316 include 'COMMON.NAMES'
1317 include 'COMMON.CHAIN'
1318 include 'COMMON.FFIELD'
1319 include 'COMMON.SBRIDGE'
1320 include 'COMMON.HEADER'
1321 include 'COMMON.CONTROL'
1322 include 'COMMON.DBASE'
1323 include 'COMMON.THREAD'
1324 include 'COMMON.TIME1'
1325 include 'COMMON.SETUP'
1326 C Read bridging residues.
1327 read (inp,*) ns,(iss(i),i=1,ns)
1329 if(me.eq.king.or..not.out1file)
1330 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1331 C Check whether the specified bridging residues are cystines.
1333 if (itype(iss(i)).ne.1) then
1334 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1335 & 'Do you REALLY think that the residue ',
1336 & restyp(itype(iss(i))),i,
1337 & ' can form a disulfide bridge?!!!'
1338 write (*,'(2a,i3,a)')
1339 & 'Do you REALLY think that the residue ',
1340 & restyp(itype(iss(i))),i,
1341 & ' can form a disulfide bridge?!!!'
1343 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1348 C Read preformed bridges.
1350 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1352 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1355 C Check if the residues involved in bridges are in the specified list of
1356 C bridging residues.
1359 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1360 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1361 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1362 & ' contains residues present in other pairs.'
1363 write (*,'(a,i3,a)') 'Disulfide pair',i,
1364 & ' contains residues present in other pairs.'
1366 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1372 if (ihpb(i).eq.iss(j)) goto 10
1374 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1377 if (jhpb(i).eq.iss(j)) goto 20
1379 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1385 ihpb(i)=ihpb(i)+nres
1386 jhpb(i)=jhpb(i)+nres
1392 c----------------------------------------------------------------------------
1393 subroutine read_x(kanal,*)
1394 implicit real*8 (a-h,o-z)
1395 include 'DIMENSIONS'
1396 include 'COMMON.GEO'
1397 include 'COMMON.VAR'
1398 include 'COMMON.CHAIN'
1399 include 'COMMON.IOUNITS'
1400 include 'COMMON.CONTROL'
1401 include 'COMMON.LOCAL'
1402 include 'COMMON.INTERACT'
1403 c Read coordinates from input
1405 read(kanal,'(8f10.5)',end=10,err=10)
1406 & ((c(l,k),l=1,3),k=1,nres),
1407 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1410 c(j,2*nres)=c(j,nres)
1412 call int_from_cart1(.false.)
1415 dc(j,i)=c(j,i+1)-c(j,i)
1416 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1420 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1422 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1423 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1431 c----------------------------------------------------------------------------
1432 subroutine read_threadbase
1433 implicit real*8 (a-h,o-z)
1434 include 'DIMENSIONS'
1435 include 'COMMON.IOUNITS'
1436 include 'COMMON.GEO'
1437 include 'COMMON.VAR'
1438 include 'COMMON.INTERACT'
1439 include 'COMMON.LOCAL'
1440 include 'COMMON.NAMES'
1441 include 'COMMON.CHAIN'
1442 include 'COMMON.FFIELD'
1443 include 'COMMON.SBRIDGE'
1444 include 'COMMON.HEADER'
1445 include 'COMMON.CONTROL'
1446 include 'COMMON.DBASE'
1447 include 'COMMON.THREAD'
1448 include 'COMMON.TIME1'
1449 C Read pattern database for threading.
1450 read (icbase,*) nseq
1452 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1453 & nres_base(2,i),nres_base(3,i)
1454 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1456 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1457 c & nres_base(2,i),nres_base(3,i)
1458 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1462 if (weidis.eq.0.0D0) weidis=0.1D0
1471 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1472 write (iout,'(a,i5)') 'nexcl: ',nexcl
1473 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1476 c------------------------------------------------------------------------------
1477 subroutine setup_var
1478 implicit real*8 (a-h,o-z)
1479 include 'DIMENSIONS'
1480 include 'COMMON.IOUNITS'
1481 include 'COMMON.GEO'
1482 include 'COMMON.VAR'
1483 include 'COMMON.INTERACT'
1484 include 'COMMON.LOCAL'
1485 include 'COMMON.NAMES'
1486 include 'COMMON.CHAIN'
1487 include 'COMMON.FFIELD'
1488 include 'COMMON.SBRIDGE'
1489 include 'COMMON.HEADER'
1490 include 'COMMON.CONTROL'
1491 include 'COMMON.DBASE'
1492 include 'COMMON.THREAD'
1493 include 'COMMON.TIME1'
1494 C Set up variable list.
1500 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1502 ialph(i,1)=nvar+nside
1506 if (indphi.gt.0) then
1508 else if (indback.gt.0) then
1513 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1516 c----------------------------------------------------------------------------
1517 subroutine gen_dist_constr
1518 C Generate CA distance constraints.
1519 implicit real*8 (a-h,o-z)
1520 include 'DIMENSIONS'
1521 include 'COMMON.IOUNITS'
1522 include 'COMMON.GEO'
1523 include 'COMMON.VAR'
1524 include 'COMMON.INTERACT'
1525 include 'COMMON.LOCAL'
1526 include 'COMMON.NAMES'
1527 include 'COMMON.CHAIN'
1528 include 'COMMON.FFIELD'
1529 include 'COMMON.SBRIDGE'
1530 include 'COMMON.HEADER'
1531 include 'COMMON.CONTROL'
1532 include 'COMMON.DBASE'
1533 include 'COMMON.THREAD'
1534 include 'COMMON.TIME1'
1535 dimension itype_pdb(maxres)
1536 common /pizda/ itype_pdb
1538 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1539 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1540 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1542 do i=nstart_sup,nstart_sup+nsup-1
1543 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1544 cd & ' seq_pdb', restyp(itype_pdb(i))
1545 do j=i+2,nstart_sup+nsup-1
1547 ihpb(nhpb)=i+nstart_seq-nstart_sup
1548 jhpb(nhpb)=j+nstart_seq-nstart_sup
1550 dhpb(nhpb)=dist(i,j)
1553 cd write (iout,'(a)') 'Distance constraints:'
1558 cd if (ii.gt.nres) then
1563 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1564 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1565 cd & dhpb(i),forcon(i)
1569 c----------------------------------------------------------------------------
1571 implicit real*8 (a-h,o-z)
1572 include 'DIMENSIONS'
1573 include 'COMMON.MAP'
1574 include 'COMMON.IOUNITS'
1575 character*3 angid(4) /'THE','PHI','ALP','OME'/
1576 character*80 mapcard,ucase
1578 read (inp,'(a)') mapcard
1579 mapcard=ucase(mapcard)
1580 if (index(mapcard,'PHI').gt.0) then
1582 else if (index(mapcard,'THE').gt.0) then
1584 else if (index(mapcard,'ALP').gt.0) then
1586 else if (index(mapcard,'OME').gt.0) then
1589 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1590 stop 'Error - illegal variable spec in MAP card.'
1592 call readi (mapcard,'RES1',res1(imap),0)
1593 call readi (mapcard,'RES2',res2(imap),0)
1594 if (res1(imap).eq.0) then
1595 res1(imap)=res2(imap)
1596 else if (res2(imap).eq.0) then
1597 res2(imap)=res1(imap)
1599 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1601 & 'Error - illegal definition of variable group in MAP.'
1602 stop 'Error - illegal definition of variable group in MAP.'
1604 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1605 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1606 call readi(mapcard,'NSTEP',nstep(imap),0)
1607 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1609 & 'Illegal boundary and/or step size specification in MAP.'
1610 stop 'Illegal boundary and/or step size specification in MAP.'
1615 c----------------------------------------------------------------------------
1617 implicit real*8 (a-h,o-z)
1618 include 'DIMENSIONS'
1619 include 'COMMON.IOUNITS'
1620 include 'COMMON.GEO'
1621 include 'COMMON.CSA'
1622 include 'COMMON.BANK'
1623 include 'COMMON.CONTROL'
1625 character*620 mcmcard
1626 call card_concat(mcmcard)
1628 call readi(mcmcard,'NCONF',nconf,50)
1629 call readi(mcmcard,'NADD',nadd,0)
1630 call readi(mcmcard,'JSTART',jstart,1)
1631 call readi(mcmcard,'JEND',jend,1)
1632 call readi(mcmcard,'NSTMAX',nstmax,500000)
1633 call readi(mcmcard,'N0',n0,1)
1634 call readi(mcmcard,'N1',n1,6)
1635 call readi(mcmcard,'N2',n2,4)
1636 call readi(mcmcard,'N3',n3,0)
1637 call readi(mcmcard,'N4',n4,0)
1638 call readi(mcmcard,'N5',n5,0)
1639 call readi(mcmcard,'N6',n6,10)
1640 call readi(mcmcard,'N7',n7,0)
1641 call readi(mcmcard,'N8',n8,0)
1642 call readi(mcmcard,'N9',n9,0)
1643 call readi(mcmcard,'N14',n14,0)
1644 call readi(mcmcard,'N15',n15,0)
1645 call readi(mcmcard,'N16',n16,0)
1646 call readi(mcmcard,'N17',n17,0)
1647 call readi(mcmcard,'N18',n18,0)
1649 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1651 call readi(mcmcard,'NDIFF',ndiff,2)
1652 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1653 call readi(mcmcard,'IS1',is1,1)
1654 call readi(mcmcard,'IS2',is2,8)
1655 call readi(mcmcard,'NRAN0',nran0,4)
1656 call readi(mcmcard,'NRAN1',nran1,2)
1657 call readi(mcmcard,'IRR',irr,1)
1658 call readi(mcmcard,'NSEED',nseed,20)
1659 call readi(mcmcard,'NTOTAL',ntotal,10000)
1660 call reada(mcmcard,'CUT1',cut1,2.0d0)
1661 call reada(mcmcard,'CUT2',cut2,5.0d0)
1662 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1663 call readi(mcmcard,'ICMAX',icmax,3)
1664 call readi(mcmcard,'IRESTART',irestart,0)
1665 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1668 call reada(mcmcard,'DELE',dele,20.0d0)
1669 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1670 call readi(mcmcard,'IREF',iref,0)
1671 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1672 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1673 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1674 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1675 write (iout,*) "NCONF_IN",nconf_in
1678 c----------------------------------------------------------------------------
1679 cfmc subroutine mcmfread
1680 cfmc implicit real*8 (a-h,o-z)
1681 cfmc include 'DIMENSIONS'
1682 cfmc include 'COMMON.MCMF'
1683 cfmc include 'COMMON.IOUNITS'
1684 cfmc include 'COMMON.GEO'
1685 cfmc character*80 ucase
1686 cfmc character*620 mcmcard
1687 cfmc call card_concat(mcmcard)
1689 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1690 cfmc write(iout,*)'MAXRANT=',maxrant
1691 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1692 cfmc write(iout,*)'MAXFAM=',maxfam
1693 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1694 cfmc write(iout,*)'NNET1=',nnet1
1695 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1696 cfmc write(iout,*)'NNET2=',nnet2
1697 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1698 cfmc write(iout,*)'NNET3=',nnet3
1699 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1700 cfmc write(iout,*)'ILASTT=',ilastt
1701 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1702 cfmc write(iout,*)'MAXSTR=',maxstr
1703 cfmc maxstr_f=maxstr/maxfam
1704 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1705 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1706 cfmc write(iout,*)'NMCMF=',nmcmf
1707 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1708 cfmc write(iout,*)'IFOCUS=',ifocus
1709 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1710 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1711 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1712 cfmc write(iout,*)'INTPRT=',intprt
1713 cfmc call readi(mcmcard,'IPRT',iprt,100)
1714 cfmc write(iout,*)'IPRT=',iprt
1715 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1716 cfmc write(iout,*)'IMAXTR=',imaxtr
1717 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1718 cfmc write(iout,*)'MAXEVEN=',maxeven
1719 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1720 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1721 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1722 cfmc write(iout,*)'INIMIN=',inimin
1723 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1724 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1725 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1726 cfmc write(iout,*)'NTHREAD=',nthread
1727 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1728 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1729 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1730 cfmc write(iout,*)'MAXPERT=',maxpert
1731 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1732 cfmc write(iout,*)'IRMSD=',irmsd
1733 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1734 cfmc write(iout,*)'DENEMIN=',denemin
1735 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1736 cfmc write(iout,*)'RCUT1S=',rcut1s
1737 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1738 cfmc write(iout,*)'RCUT1E=',rcut1e
1739 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1740 cfmc write(iout,*)'RCUT2S=',rcut2s
1741 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1742 cfmc write(iout,*)'RCUT2E=',rcut2e
1743 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1744 cfmc write(iout,*)'DPERT1=',d_pert1
1745 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1746 cfmc write(iout,*)'DPERT1A=',d_pert1a
1747 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1748 cfmc write(iout,*)'DPERT2=',d_pert2
1749 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1750 cfmc write(iout,*)'DPERT2A=',d_pert2a
1751 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1752 cfmc write(iout,*)'DPERT2B=',d_pert2b
1753 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1754 cfmc write(iout,*)'DPERT2C=',d_pert2c
1755 cfmc d_pert1=deg2rad*d_pert1
1756 cfmc d_pert1a=deg2rad*d_pert1a
1757 cfmc d_pert2=deg2rad*d_pert2
1758 cfmc d_pert2a=deg2rad*d_pert2a
1759 cfmc d_pert2b=deg2rad*d_pert2b
1760 cfmc d_pert2c=deg2rad*d_pert2c
1761 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1762 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1763 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1764 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1765 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1766 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1767 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1768 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1769 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1770 cfmc write(iout,*)'RCUTINI=',rcutini
1771 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1772 cfmc write(iout,*)'GRAT=',grat
1773 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1774 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1778 c----------------------------------------------------------------------------
1780 implicit real*8 (a-h,o-z)
1781 include 'DIMENSIONS'
1782 include 'COMMON.MCM'
1783 include 'COMMON.MCE'
1784 include 'COMMON.IOUNITS'
1786 character*320 mcmcard
1787 call card_concat(mcmcard)
1788 call readi(mcmcard,'MAXACC',maxacc,100)
1789 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1790 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1791 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1792 call readi(mcmcard,'MAXREPM',maxrepm,200)
1793 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1794 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1795 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1796 call reada(mcmcard,'E_UP',e_up,5.0D0)
1797 call reada(mcmcard,'DELTE',delte,0.1D0)
1798 call readi(mcmcard,'NSWEEP',nsweep,5)
1799 call readi(mcmcard,'NSTEPH',nsteph,0)
1800 call readi(mcmcard,'NSTEPC',nstepc,0)
1801 call reada(mcmcard,'TMIN',tmin,298.0D0)
1802 call reada(mcmcard,'TMAX',tmax,298.0D0)
1803 call readi(mcmcard,'NWINDOW',nwindow,0)
1804 call readi(mcmcard,'PRINT_MC',print_mc,0)
1805 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1806 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1807 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1808 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1809 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1810 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1811 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1812 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1813 if (nwindow.gt.0) then
1814 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1816 winlen(i)=winend(i)-winstart(i)+1
1819 if (tmax.lt.tmin) tmax=tmin
1820 if (tmax.eq.tmin) then
1824 if (nstepc.gt.0 .and. nsteph.gt.0) then
1825 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1826 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1828 C Probabilities of different move types
1829 sumpro_type(0)=0.0D0
1830 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1831 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1832 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1833 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1834 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1835 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1836 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1838 print *,'i',i,' sumprotype',sumpro_type(i)
1839 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1840 print *,'i',i,' sumprotype',sumpro_type(i)
1844 c----------------------------------------------------------------------------
1845 subroutine read_minim
1846 implicit real*8 (a-h,o-z)
1847 include 'DIMENSIONS'
1848 include 'COMMON.MINIM'
1849 include 'COMMON.IOUNITS'
1850 include 'COMMON.CONTROL'
1851 include 'COMMON.SETUP'
1853 character*320 minimcard
1854 call card_concat(minimcard)
1855 call readi(minimcard,'MAXMIN',maxmin,2000)
1856 call readi(minimcard,'MAXFUN',maxfun,5000)
1857 call readi(minimcard,'MINMIN',minmin,maxmin)
1858 call readi(minimcard,'MINFUN',minfun,maxmin)
1859 call reada(minimcard,'TOLF',tolf,1.0D-2)
1860 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1861 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1862 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1863 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1865 if (.not. out1file .or. me.eq.king) then
1867 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1868 & 'Options in energy minimization:'
1869 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1870 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1871 & 'MinMin:',MinMin,' MinFun:',MinFun,
1872 & ' TolF:',TolF,' RTolF:',RTolF
1878 c----------------------------------------------------------------------------
1879 subroutine read_angles(kanal,*)
1880 implicit real*8 (a-h,o-z)
1881 include 'DIMENSIONS'
1882 include 'COMMON.GEO'
1883 include 'COMMON.VAR'
1884 include 'COMMON.CHAIN'
1885 include 'COMMON.IOUNITS'
1886 include 'COMMON.CONTROL'
1887 c Read angles from input
1889 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1890 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1891 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1892 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1895 c 9/7/01 avoid 180 deg valence angle
1896 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1898 theta(i)=deg2rad*theta(i)
1899 phi(i)=deg2rad*phi(i)
1900 alph(i)=deg2rad*alph(i)
1901 omeg(i)=deg2rad*omeg(i)
1906 c----------------------------------------------------------------------------
1907 subroutine reada(rekord,lancuch,wartosc,default)
1909 character*(*) rekord,lancuch
1910 double precision wartosc,default
1913 iread=index(rekord,lancuch)
1914 if (iread.eq.0) then
1918 iread=iread+ilen(lancuch)+1
1919 read (rekord(iread:),*,err=10,end=10) wartosc
1924 c----------------------------------------------------------------------------
1925 subroutine readi(rekord,lancuch,wartosc,default)
1927 character*(*) rekord,lancuch
1928 integer wartosc,default
1931 iread=index(rekord,lancuch)
1932 if (iread.eq.0) then
1936 iread=iread+ilen(lancuch)+1
1937 read (rekord(iread:),*,err=10,end=10) wartosc
1942 c----------------------------------------------------------------------------
1943 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1946 integer tablica(dim),default
1947 character*(*) rekord,lancuch
1954 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1955 if (iread.eq.0) return
1956 iread=iread+ilen(lancuch)+1
1957 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1960 c----------------------------------------------------------------------------
1961 subroutine multreada(rekord,lancuch,tablica,dim,default)
1964 double precision tablica(dim),default
1965 character*(*) rekord,lancuch
1972 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1973 if (iread.eq.0) return
1974 iread=iread+ilen(lancuch)+1
1975 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1978 c----------------------------------------------------------------------------
1979 subroutine openunits
1980 implicit real*8 (a-h,o-z)
1981 include 'DIMENSIONS'
1984 character*16 form,nodename
1987 include 'COMMON.SETUP'
1988 include 'COMMON.IOUNITS'
1990 include 'COMMON.CONTROL'
1991 integer lenpre,lenpot,ilen,lentmp
1993 character*3 out1file_text,ucase
1996 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1997 call getenv_loc("PREFIX",prefix)
1999 call getenv_loc("POT",pot)
2000 call getenv_loc("DIRTMP",tmpdir)
2001 call getenv_loc("CURDIR",curdir)
2002 call getenv_loc("OUT1FILE",out1file_text)
2003 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2004 out1file_text=ucase(out1file_text)
2005 if (out1file_text(1:1).eq."Y") then
2008 out1file=fg_rank.gt.0
2013 if (lentmp.gt.0) then
2014 write (*,'(80(1h!))')
2015 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2016 write (*,'(80(1h!))')
2017 write (*,*)"All output files will be on node /tmp directory."
2019 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2020 if (me.eq.king) then
2021 write (*,*) "The master node is ",nodename
2022 else if (fg_rank.eq.0) then
2023 write (*,*) "I am the CG slave node ",nodename
2025 write (*,*) "I am the FG slave node ",nodename
2028 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2029 lenpre = lentmp+lenpre+1
2031 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2032 C Get the names and open the input files
2033 #if defined(WINIFL) || defined(WINPGI)
2034 open(1,file=pref_orig(:ilen(pref_orig))//
2035 & '.inp',status='old',readonly,shared)
2036 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2037 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2038 C Get parameter filenames and open the parameter files.
2039 call getenv_loc('BONDPAR',bondname)
2040 open (ibond,file=bondname,status='old',readonly,shared)
2041 call getenv_loc('THETPAR',thetname)
2042 open (ithep,file=thetname,status='old',readonly,shared)
2044 call getenv_loc('THETPARPDB',thetname_pdb)
2045 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2047 call getenv_loc('ROTPAR',rotname)
2048 open (irotam,file=rotname,status='old',readonly,shared)
2050 call getenv_loc('ROTPARPDB',rotname_pdb)
2051 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2053 call getenv_loc('TORPAR',torname)
2054 open (itorp,file=torname,status='old',readonly,shared)
2055 call getenv_loc('TORDPAR',tordname)
2056 open (itordp,file=tordname,status='old',readonly,shared)
2057 call getenv_loc('FOURIER',fouriername)
2058 open (ifourier,file=fouriername,status='old',readonly,shared)
2059 call getenv_loc('ELEPAR',elename)
2060 open (ielep,file=elename,status='old',readonly,shared)
2061 call getenv_loc('SIDEPAR',sidename)
2062 open (isidep,file=sidename,status='old',readonly,shared)
2063 #elif (defined CRAY) || (defined AIX)
2064 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2066 c print *,"Processor",myrank," opened file 1"
2067 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2068 c print *,"Processor",myrank," opened file 9"
2069 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2070 C Get parameter filenames and open the parameter files.
2071 call getenv_loc('BONDPAR',bondname)
2072 open (ibond,file=bondname,status='old',action='read')
2073 c print *,"Processor",myrank," opened file IBOND"
2074 call getenv_loc('THETPAR',thetname)
2075 open (ithep,file=thetname,status='old',action='read')
2076 c print *,"Processor",myrank," opened file ITHEP"
2078 call getenv_loc('THETPARPDB',thetname_pdb)
2079 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2081 call getenv_loc('ROTPAR',rotname)
2082 open (irotam,file=rotname,status='old',action='read')
2083 c print *,"Processor",myrank," opened file IROTAM"
2085 call getenv_loc('ROTPARPDB',rotname_pdb)
2086 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2088 call getenv_loc('TORPAR',torname)
2089 open (itorp,file=torname,status='old',action='read')
2090 c print *,"Processor",myrank," opened file ITORP"
2091 call getenv_loc('TORDPAR',tordname)
2092 open (itordp,file=tordname,status='old',action='read')
2093 c print *,"Processor",myrank," opened file ITORDP"
2094 call getenv_loc('SCCORPAR',sccorname)
2095 open (isccor,file=sccorname,status='old',action='read')
2096 c print *,"Processor",myrank," opened file ISCCOR"
2097 call getenv_loc('FOURIER',fouriername)
2098 open (ifourier,file=fouriername,status='old',action='read')
2099 c print *,"Processor",myrank," opened file IFOURIER"
2100 call getenv_loc('ELEPAR',elename)
2101 open (ielep,file=elename,status='old',action='read')
2102 c print *,"Processor",myrank," opened file IELEP"
2103 call getenv_loc('SIDEPAR',sidename)
2104 open (isidep,file=sidename,status='old',action='read')
2105 c print *,"Processor",myrank," opened file ISIDEP"
2106 c print *,"Processor",myrank," opened parameter files"
2108 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2109 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2110 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2111 C Get parameter filenames and open the parameter files.
2112 call getenv_loc('BONDPAR',bondname)
2113 open (ibond,file=bondname,status='old')
2114 call getenv_loc('THETPAR',thetname)
2115 open (ithep,file=thetname,status='old')
2117 call getenv_loc('THETPARPDB',thetname_pdb)
2118 open (ithep_pdb,file=thetname_pdb,status='old')
2120 call getenv_loc('ROTPAR',rotname)
2121 open (irotam,file=rotname,status='old')
2123 call getenv_loc('ROTPARPDB',rotname_pdb)
2124 open (irotam_pdb,file=rotname_pdb,status='old')
2126 call getenv_loc('TORPAR',torname)
2127 open (itorp,file=torname,status='old')
2128 call getenv_loc('TORDPAR',tordname)
2129 open (itordp,file=tordname,status='old')
2130 call getenv_loc('SCCORPAR',sccorname)
2131 open (isccor,file=sccorname,status='old')
2132 call getenv_loc('FOURIER',fouriername)
2133 open (ifourier,file=fouriername,status='old')
2134 call getenv_loc('ELEPAR',elename)
2135 open (ielep,file=elename,status='old')
2136 call getenv_loc('SIDEPAR',sidename)
2137 open (isidep,file=sidename,status='old')
2138 call getenv_loc('LIPTRANPAR',liptranname)
2139 open (iliptranpar,file=liptranname,status='old')
2141 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2143 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2144 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2145 C Get parameter filenames and open the parameter files.
2146 call getenv_loc('BONDPAR',bondname)
2147 open (ibond,file=bondname,status='old',readonly)
2148 call getenv_loc('THETPAR',thetname)
2149 open (ithep,file=thetname,status='old',readonly)
2150 call getenv_loc('ROTPAR',rotname)
2151 open (irotam,file=rotname,status='old',readonly)
2152 call getenv_loc('TORPAR',torname)
2153 open (itorp,file=torname,status='old',readonly)
2154 call getenv_loc('TORDPAR',tordname)
2155 open (itordp,file=tordname,status='old',readonly)
2156 call getenv_loc('SCCORPAR',sccorname)
2157 open (isccor,file=sccorname,status='old',readonly)
2159 call getenv_loc('THETPARPDB',thetname_pdb)
2160 c print *,"thetname_pdb ",thetname_pdb
2161 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2162 c print *,ithep_pdb," opened"
2164 call getenv_loc('FOURIER',fouriername)
2165 open (ifourier,file=fouriername,status='old',readonly)
2166 call getenv_loc('ELEPAR',elename)
2167 open (ielep,file=elename,status='old',readonly)
2168 call getenv_loc('SIDEPAR',sidename)
2169 open (isidep,file=sidename,status='old',readonly)
2170 call getenv_loc('LIPTRANPAR',liptranname)
2171 open (iliptranpar,file=liptranname,status='old',readonly)
2173 call getenv_loc('ROTPARPDB',rotname_pdb)
2174 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2179 C 8/9/01 In the newest version SCp interaction constants are read from a file
2180 C Use -DOLDSCP to use hard-coded constants instead.
2182 call getenv_loc('SCPPAR',scpname)
2183 #if defined(WINIFL) || defined(WINPGI)
2184 open (iscpp,file=scpname,status='old',readonly,shared)
2185 #elif (defined CRAY) || (defined AIX)
2186 open (iscpp,file=scpname,status='old',action='read')
2188 open (iscpp,file=scpname,status='old')
2190 open (iscpp,file=scpname,status='old',readonly)
2193 call getenv_loc('PATTERN',patname)
2194 #if defined(WINIFL) || defined(WINPGI)
2195 open (icbase,file=patname,status='old',readonly,shared)
2196 #elif (defined CRAY) || (defined AIX)
2197 open (icbase,file=patname,status='old',action='read')
2199 open (icbase,file=patname,status='old')
2201 open (icbase,file=patname,status='old',readonly)
2204 C Open output file only for CG processes
2205 c print *,"Processor",myrank," fg_rank",fg_rank
2206 if (fg_rank.eq.0) then
2208 if (nodes.eq.1) then
2211 npos = dlog10(dfloat(nodes-1))+1
2213 if (npos.lt.3) npos=3
2214 write (liczba,'(i1)') npos
2215 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2217 write (liczba,form) me
2218 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2219 & liczba(:ilen(liczba))
2220 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2222 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2224 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2225 & liczba(:ilen(liczba))//'.mol2'
2226 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2227 & liczba(:ilen(liczba))//'.stat'
2229 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2230 & //liczba(:ilen(liczba))//'.stat')
2231 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2234 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2235 & liczba(:ilen(liczba))//'.const'
2240 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2241 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2242 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2243 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2244 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2246 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2248 rest2name=prefix(:ilen(prefix))//'.rst'
2250 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2253 #if defined(AIX) || defined(PGI)
2254 if (me.eq.king .or. .not. out1file)
2255 & open(iout,file=outname,status='unknown')
2257 if (fg_rank.gt.0) then
2258 write (liczba,'(i3.3)') myrank/nfgtasks
2259 write (ll,'(bz,i3.3)') fg_rank
2260 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2265 open(igeom,file=intname,status='unknown',position='append')
2266 open(ipdb,file=pdbname,status='unknown')
2267 open(imol2,file=mol2name,status='unknown')
2268 open(istat,file=statname,status='unknown',position='append')
2270 c1out open(iout,file=outname,status='unknown')
2273 if (me.eq.king .or. .not.out1file)
2274 & open(iout,file=outname,status='unknown')
2276 if (fg_rank.gt.0) then
2277 write (liczba,'(i3.3)') myrank/nfgtasks
2278 write (ll,'(bz,i3.3)') fg_rank
2279 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2284 open(igeom,file=intname,status='unknown',access='append')
2285 open(ipdb,file=pdbname,status='unknown')
2286 open(imol2,file=mol2name,status='unknown')
2287 open(istat,file=statname,status='unknown',access='append')
2289 c1out open(iout,file=outname,status='unknown')
2292 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2293 csa_seed=prefix(:lenpre)//'.CSA.seed'
2294 csa_history=prefix(:lenpre)//'.CSA.history'
2295 csa_bank=prefix(:lenpre)//'.CSA.bank'
2296 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2297 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2298 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2299 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2300 csa_int=prefix(:lenpre)//'.int'
2301 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2302 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2303 csa_in=prefix(:lenpre)//'.CSA.in'
2304 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2307 write (iout,'(80(1h-))')
2308 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2309 write (iout,'(80(1h-))')
2310 write (iout,*) "Input file : ",
2311 & pref_orig(:ilen(pref_orig))//'.inp'
2312 write (iout,*) "Output file : ",
2313 & outname(:ilen(outname))
2315 write (iout,*) "Sidechain potential file : ",
2316 & sidename(:ilen(sidename))
2318 write (iout,*) "SCp potential file : ",
2319 & scpname(:ilen(scpname))
2321 write (iout,*) "Electrostatic potential file : ",
2322 & elename(:ilen(elename))
2323 write (iout,*) "Cumulant coefficient file : ",
2324 & fouriername(:ilen(fouriername))
2325 write (iout,*) "Torsional parameter file : ",
2326 & torname(:ilen(torname))
2327 write (iout,*) "Double torsional parameter file : ",
2328 & tordname(:ilen(tordname))
2329 write (iout,*) "SCCOR parameter file : ",
2330 & sccorname(:ilen(sccorname))
2331 write (iout,*) "Bond & inertia constant file : ",
2332 & bondname(:ilen(bondname))
2333 write (iout,*) "Bending parameter file : ",
2334 & thetname(:ilen(thetname))
2335 write (iout,*) "Rotamer parameter file : ",
2336 & rotname(:ilen(rotname))
2337 write (iout,*) "Thetpdb parameter file : ",
2338 & thetname_pdb(:ilen(thetname_pdb))
2339 write (iout,*) "Threading database : ",
2340 & patname(:ilen(patname))
2342 &write (iout,*)" DIRTMP : ",
2344 write (iout,'(80(1h-))')
2348 c----------------------------------------------------------------------------
2349 subroutine card_concat(card)
2350 implicit real*8 (a-h,o-z)
2351 include 'DIMENSIONS'
2352 include 'COMMON.IOUNITS'
2354 character*80 karta,ucase
2356 read (inp,'(a)') karta
2359 do while (karta(80:80).eq.'&')
2360 card=card(:ilen(card)+1)//karta(:79)
2361 read (inp,'(a)') karta
2364 card=card(:ilen(card)+1)//karta
2367 c----------------------------------------------------------------------------------
2369 implicit real*8 (a-h,o-z)
2370 include 'DIMENSIONS'
2371 include 'COMMON.CHAIN'
2372 include 'COMMON.IOUNITS'
2374 include 'COMMON.CONTROL'
2375 open(irest2,file=rest2name,status='unknown')
2376 read(irest2,*) totT,EK,potE,totE,t_bath
2379 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2382 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2384 if(usampl.or.homol_nset.gt.1) then
2385 read (irest2,*) iset
2390 c---------------------------------------------------------------------------------
2391 subroutine read_fragments
2392 implicit real*8 (a-h,o-z)
2393 include 'DIMENSIONS'
2397 include 'COMMON.SETUP'
2398 include 'COMMON.CHAIN'
2399 include 'COMMON.IOUNITS'
2401 include 'COMMON.CONTROL'
2403 read(inp,*) nset,nfrag,npair,nfrag_back
2404 if(me.eq.king.or..not.out1file)
2405 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2406 & " nfrag_back",nfrag_back
2408 read(inp,*) mset(iset1)
2410 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2412 if(me.eq.king.or..not.out1file)
2413 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2414 & ifrag(2,i,iset1), qinfrag(i,iset1)
2417 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2419 if(me.eq.king.or..not.out1file)
2420 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2421 & ipair(2,i,iset1), qinpair(i,iset1)
2424 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2425 & wfrag_back(3,i,iset1),
2426 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2427 if(me.eq.king.or..not.out1file)
2428 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2429 & wfrag_back(2,i,iset1),
2430 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2431 & ifrag_back(2,i,iset1)
2436 C---------------------------------------------------------------------------
2437 subroutine read_afminp
2438 implicit real*8 (a-h,o-z)
2439 include 'DIMENSIONS'
2443 include 'COMMON.SETUP'
2444 include 'COMMON.CONTROL'
2445 include 'COMMON.CHAIN'
2446 include 'COMMON.IOUNITS'
2447 include 'COMMON.SBRIDGE'
2448 character*320 afmcard
2449 c print *, "wchodze"
2450 call card_concat(afmcard)
2451 call readi(afmcard,"BEG",afmbeg,0)
2452 call readi(afmcard,"END",afmend,0)
2453 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2454 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2455 print *,'FORCE=' ,forceAFMconst
2456 CCCC NOW PROPERTIES FOR AFM
2459 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2461 distafminit=dsqrt(distafminit)
2462 print *,'initdist',distafminit
2466 c-------------------------------------------------------------------------------
2467 subroutine read_saxs_constr
2468 implicit real*8 (a-h,o-z)
2469 include 'DIMENSIONS'
2473 include 'COMMON.SETUP'
2474 include 'COMMON.CONTROL'
2475 include 'COMMON.CHAIN'
2476 include 'COMMON.IOUNITS'
2477 include 'COMMON.SBRIDGE'
2478 double precision cm(3)
2480 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2482 if (saxs_mode.eq.0) then
2483 c SAXS distance distribution
2485 read(inp,*) distsaxs(i),Psaxs(i)
2489 Cnorm = Cnorm + Psaxs(i)
2491 write (iout,*) "Cnorm",Cnorm
2493 Psaxs(i)=Psaxs(i)/Cnorm
2495 write (iout,*) "Normalized distance distribution from SAXS"
2497 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2502 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2509 cm(j)=cm(j)+Csaxs(j,i)
2517 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2520 write (iout,*) "SAXS sphere coordinates"
2522 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2527 c-------------------------------------------------------------------------------
2528 subroutine read_dist_constr
2529 implicit real*8 (a-h,o-z)
2530 include 'DIMENSIONS'
2534 include 'COMMON.SETUP'
2535 include 'COMMON.CONTROL'
2536 include 'COMMON.CHAIN'
2537 include 'COMMON.IOUNITS'
2538 include 'COMMON.SBRIDGE'
2539 integer ifrag_(2,100),ipair_(2,100)
2540 double precision wfrag_(100),wpair_(100)
2541 character*500 controlcard
2542 c write (iout,*) "Calling read_dist_constr"
2543 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2545 call card_concat(controlcard)
2546 call readi(controlcard,"NFRAG",nfrag_,0)
2547 call readi(controlcard,"NPAIR",npair_,0)
2548 call readi(controlcard,"NDIST",ndist_,0)
2549 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2550 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2551 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2552 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2553 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2554 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2555 c write (iout,*) "IFRAG"
2557 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2559 c write (iout,*) "IPAIR"
2561 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2563 if (.not.refstr .and. nfrag.gt.0) then
2565 & "ERROR: no reference structure to compute distance restraints"
2567 & "Restraints must be specified explicitly (NDIST=number)"
2570 if (nfrag.lt.2 .and. npair.gt.0) then
2571 write (iout,*) "ERROR: Less than 2 fragments specified",
2572 & " but distance restraints between pairs requested"
2577 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2578 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2579 & ifrag_(2,i)=nstart_sup+nsup-1
2580 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2582 if (wfrag_(i).gt.0.0d0) then
2583 do j=ifrag_(1,i),ifrag_(2,i)-1
2584 do k=j+1,ifrag_(2,i)
2585 c write (iout,*) "j",j," k",k
2587 if (constr_dist.eq.1) then
2592 forcon(nhpb)=wfrag_(i)
2593 else if (constr_dist.eq.2) then
2594 if (ddjk.le.dist_cut) then
2599 forcon(nhpb)=wfrag_(i)
2606 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2609 if (.not.out1file .or. me.eq.king)
2610 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2611 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2613 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2614 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2621 if (wpair_(i).gt.0.0d0) then
2629 do j=ifrag_(1,ii),ifrag_(2,ii)
2630 do k=ifrag_(1,jj),ifrag_(2,jj)
2634 forcon(nhpb)=wpair_(i)
2635 dhpb(nhpb)=dist(j,k)
2637 if (.not.out1file .or. me.eq.king)
2638 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2639 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2641 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2642 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2649 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2650 if (forcon(nhpb+1).gt.0.0d0) then
2652 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2654 if (.not.out1file .or. me.eq.king)
2655 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2656 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2658 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2659 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2666 c-------------------------------------------------------------------------------
2668 subroutine read_constr_homology
2670 include 'DIMENSIONS'
2674 include 'COMMON.SETUP'
2675 include 'COMMON.CONTROL'
2676 include 'COMMON.CHAIN'
2677 include 'COMMON.IOUNITS'
2679 include 'COMMON.GEO'
2680 include 'COMMON.INTERACT'
2681 include 'COMMON.NAMES'
2683 c For new homol impl
2685 include 'COMMON.VAR'
2688 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2690 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2691 c & sigma_odl_temp(maxres,maxres,max_template)
2693 character*24 model_ki_dist, model_ki_angle
2694 character*500 controlcard
2695 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2696 logical lprn /.true./
2701 c FP - Nov. 2014 Temporary specifications for new vars
2703 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2705 double precision, dimension (max_template,maxres) :: rescore
2706 double precision, dimension (max_template,maxres) :: rescore2
2707 double precision, dimension (max_template,maxres) :: rescore3
2708 character*24 pdbfile,tpl_k_rescore
2709 c -----------------------------------------------------------------
2710 c Reading multiple PDB ref structures and calculation of retraints
2711 c not using pre-computed ones stored in files model_ki_{dist,angle}
2713 c -----------------------------------------------------------------
2716 c Alternative: reading from input
2717 call card_concat(controlcard)
2718 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2719 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2720 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2721 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2722 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2723 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2724 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2725 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2726 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2727 if(.not.read2sigma.and.start_from_model) then
2728 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
2729 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2730 start_from_model=.false.
2732 if(start_from_model .and. (me.eq.king .or. .not. out1file))
2733 & write(iout,*) 'START_FROM_MODELS is ON'
2734 if(start_from_model .and. rest) then
2735 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2736 write(iout,*) 'START_FROM_MODELS is OFF'
2737 write(iout,*) 'remove restart keyword from input'
2740 if (homol_nset.gt.1)then
2741 call card_concat(controlcard)
2742 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2743 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2744 write(iout,*) "iset homology_weight "
2746 write(iout,*) i,waga_homology(i)
2749 iset=mod(kolor,homol_nset)+1
2752 waga_homology(1)=1.0
2755 cd write (iout,*) "nnt",nnt," nct",nct
2762 c write(iout,*) 'nnt=',nnt,'nct=',nct
2765 do k=1,constr_homology
2778 do k=1,constr_homology
2780 read(inp,'(a)') pdbfile
2781 if(me.eq.king .or. .not. out1file)
2782 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2783 & pdbfile(:ilen(pdbfile))
2784 open(ipdbin,file=pdbfile,status='old',err=33)
2786 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
2787 & pdbfile(:ilen(pdbfile))
2790 c print *,'Begin reading pdb data'
2792 c Files containing res sim or local scores (former containing sigmas)
2795 write(kic2,'(bz,i2.2)') k
2797 tpl_k_rescore="template"//kic2//".sco"
2800 if (read2sigma) then
2801 call readpdb_template(k)
2806 c Distance restraints
2809 C Copy the coordinates from reference coordinates (?)
2813 c write (iout,*) "c(",j,i,") =",c(j,i)
2817 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2819 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2820 open (ientin,file=tpl_k_rescore,status='old')
2821 if (nnt.gt.1) rescore(k,1)=0.0d0
2822 do irec=nnt,nct ! loop for reading res sim
2823 if (read2sigma) then
2824 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
2825 & rescore3_tmp,idomain_tmp
2827 idomain(k,i_tmp)=idomain_tmp
2828 rescore(k,i_tmp)=rescore_tmp
2829 rescore2(k,i_tmp)=rescore2_tmp
2830 rescore3(k,i_tmp)=rescore3_tmp
2831 if (.not. out1file .or. me.eq.king)
2832 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
2833 & i_tmp,rescore2_tmp,rescore_tmp,
2834 & rescore3_tmp,idomain_tmp
2837 read (ientin,*,end=1401) rescore_tmp
2839 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2840 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2841 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2846 if (waga_dist.ne.0.0d0) then
2854 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2855 c write (iout,*) k,i,j,distal,dist2_cut
2857 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2858 & .and. distal.le.dist2_cut ) then
2864 c write (iout,*) "k",k
2865 c write (iout,*) "i",i," j",j," constr_homology",
2870 if (read2sigma) then
2873 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2875 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2876 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
2877 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2879 if (odl(k,ii).le.dist_cut) then
2880 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2883 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2884 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2886 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2887 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2891 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2894 l_homo(k,ii)=.false.
2901 c Theta, dihedral and SC retraints
2903 if (waga_angle.gt.0.0d0) then
2904 c open (ientin,file=tpl_k_sigma_dih,status='old')
2905 c do irec=1,maxres-3 ! loop for reading sigma_dih
2906 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2907 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2908 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2909 c & sigma_dih(k,i+nnt-1)
2914 if (idomain(k,i).eq.0) then
2918 dih(k,i)=phiref(i) ! right?
2919 c read (ientin,*) sigma_dih(k,i) ! original variant
2920 c write (iout,*) "dih(",k,i,") =",dih(k,i)
2921 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2922 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2923 c & "rescore(",k,i-2,") =",rescore(k,i-2),
2924 c & "rescore(",k,i-3,") =",rescore(k,i-3)
2926 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
2927 & rescore(k,i-2)+rescore(k,i-3))/4.0
2928 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2929 c write (iout,*) "Raw sigmas for dihedral angle restraints"
2930 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2931 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2932 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
2933 c Instead of res sim other local measure of b/b str reliability possible
2934 if (sigma_dih(k,i).ne.0)
2935 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2936 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2941 if (waga_theta.gt.0.0d0) then
2942 c open (ientin,file=tpl_k_sigma_theta,status='old')
2943 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2944 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2945 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2946 c & sigma_theta(k,i+nnt-1)
2951 do i = nnt+2,nct ! right? without parallel.
2952 c do i = i=1,nres ! alternative for bounds acc to readpdb?
2953 c do i=ithet_start,ithet_end ! with FG parallel.
2954 if (idomain(k,i).eq.0) then
2955 sigma_theta(k,i)=0.0
2958 thetatpl(k,i)=thetaref(i)
2959 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2960 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2961 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2962 c & "rescore(",k,i-2,") =",rescore(k,i-2)
2963 c read (ientin,*) sigma_theta(k,i) ! 1st variant
2964 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
2965 & rescore(k,i-2))/3.0
2966 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
2967 if (sigma_theta(k,i).ne.0)
2968 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2970 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2971 c rescore(k,i-2) ! right expression ?
2972 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2976 if (waga_d.gt.0.0d0) then
2977 c open (ientin,file=tpl_k_sigma_d,status='old')
2978 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2979 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2980 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2981 c & sigma_d(k,i+nnt-1)
2985 do i = nnt,nct ! right? without parallel.
2986 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
2987 c do i=loc_start,loc_end ! with FG parallel.
2988 if (itype(i).eq.10) cycle
2989 if (idomain(k,i).eq.0 ) then
2996 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2997 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2998 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2999 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3000 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3001 if (sigma_d(k,i).ne.0)
3002 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3004 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3005 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3006 c read (ientin,*) sigma_d(k,i) ! 1st variant
3011 c remove distance restraints not used in any model from the list
3012 c shift data in all arrays
3014 if (waga_dist.ne.0.0d0) then
3020 if (ii_in_use(ii).eq.0.and.liiflag) then
3024 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3025 & .not.liiflag.and.ii.eq.lim_odl) then
3026 if (ii.eq.lim_odl) then
3027 iishift=ii-iistart+1
3032 do ki=iistart,lim_odl-iishift
3033 ires_homo(ki)=ires_homo(ki+iishift)
3034 jres_homo(ki)=jres_homo(ki+iishift)
3035 ii_in_use(ki)=ii_in_use(ki+iishift)
3036 do k=1,constr_homology
3037 odl(k,ki)=odl(k,ki+iishift)
3038 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3039 l_homo(k,ki)=l_homo(k,ki+iishift)
3043 lim_odl=lim_odl-iishift
3048 if (constr_homology.gt.0) call homology_partition
3049 if (constr_homology.gt.0) call init_int_table
3050 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3051 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3055 if (.not.lprn) return
3056 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3057 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3058 write (iout,*) "Distance restraints from templates"
3060 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3061 & ii,ires_homo(ii),jres_homo(ii),
3062 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3063 & ki=1,constr_homology)
3065 write (iout,*) "Dihedral angle restraints from templates"
3067 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3068 & (rad2deg*dih(ki,i),
3069 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3071 write (iout,*) "Virtual-bond angle restraints from templates"
3073 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3074 & (rad2deg*thetatpl(ki,i),
3075 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3077 write (iout,*) "SC restraints from templates"
3079 write(iout,'(i5,100(4f8.2,4x))') i,
3080 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3081 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3084 c -----------------------------------------------------------------
3087 c----------------------------------------------------------------------
3090 subroutine flush(iu)
3095 subroutine flush(iu)
3100 c------------------------------------------------------------------------------
3101 subroutine copy_to_tmp(source)
3102 include "DIMENSIONS"
3103 include "COMMON.IOUNITS"
3104 character*(*) source
3105 character* 256 tmpfile
3109 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3110 inquire(file=tmpfile,exist=ex)
3112 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3113 & " to temporary directory..."
3114 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3115 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3119 c------------------------------------------------------------------------------
3120 subroutine move_from_tmp(source)
3121 include "DIMENSIONS"
3122 include "COMMON.IOUNITS"
3123 character*(*) source
3126 write (*,*) "Moving ",source(:ilen(source)),
3127 & " from temporary directory to working directory"
3128 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3129 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3132 c------------------------------------------------------------------------------
3133 subroutine random_init(seed)
3135 C Initialize random number generator
3137 implicit real*8 (a-h,o-z)
3138 include 'DIMENSIONS'
3144 logical OKRandom, prng_restart
3146 integer iseed_array(4)
3148 include 'COMMON.IOUNITS'
3149 include 'COMMON.TIME1'
3150 include 'COMMON.THREAD'
3151 include 'COMMON.SBRIDGE'
3152 include 'COMMON.CONTROL'
3153 include 'COMMON.MCM'
3154 include 'COMMON.MAP'
3155 include 'COMMON.HEADER'
3156 include 'COMMON.CSA'
3157 include 'COMMON.CHAIN'
3158 include 'COMMON.MUCA'
3160 include 'COMMON.FFIELD'
3161 include 'COMMON.SETUP'
3162 iseed=-dint(dabs(seed))
3163 if (iseed.eq.0) then
3164 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3165 & 'Random seed undefined. The program will stop.'
3166 write (*,'(/80(1h*)/20x,a/80(1h*))')
3167 & 'Random seed undefined. The program will stop.'
3169 call mpi_finalize(mpi_comm_world,ierr)
3171 stop 'Bad random seed.'
3174 if (fg_rank.eq.0) then
3178 if(me.eq.king .or. .not. out1file)
3179 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3180 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3181 OKRandom = prng_restart(me,iseedi8)
3184 tmp=65536.0d0**(4-i)
3185 iseed_array(i) = dint(seed/tmp)
3186 seed=seed-iseed_array(i)*tmp
3189 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3190 & (iseed_array(i),i=1,4)
3191 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3192 & (iseed_array(i),i=1,4)
3193 OKRandom = prng_restart(me,iseed_array)
3196 c r1 = prng_next(me)
3197 r1=ran_number(0.0D0,1.0D0)
3198 if(me.eq.king .or. .not. out1file)
3199 & write (iout,*) 'ran_num',r1
3200 if (r1.lt.0.0d0) OKRandom=.false.
3202 if (.not.OKRandom) then
3203 write (iout,*) 'PRNG IS NOT WORKING!!!'
3204 print *,'PRNG IS NOT WORKING!!!'
3207 call mpi_abort(mpi_comm_world,error_msg,ierr)
3210 write (iout,*) 'too many processors for parallel prng'
3211 write (*,*) 'too many processors for parallel prng'
3219 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)